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Abstract 


Hierarchically structured active materials in lithium-ion battery (LIB) electrodes 
with porous secondary particles are promising candidates to increase gravimetric 
energy density and rate performance of the cell. However, there are still aspects 
to this technology which are not fully understood. In order to obtain deeper 
knowledge, the goal of this work is to develop efficient tools to compute effective 
transport properties of granular cathode structures as well as porous secondary 
particles and to evaluate their influence on the electrochemical performance of 
the whole cell. 


On the one hand, the resistor network method (RN)—a tool to efficiently 
compute the effective transport properties of particulate systems—is extended 
with regard to the transport through the solid and the pore phase of granular 
media represented by sphere packings with polydisperse size-distributions. As 
for the solid phase, transport through the volume, via the surface, or a mix of 
both are considered. For all cases, appropriate analytically derived formulas 
from literature describing resistance between two single particles are used or 
combined accordingly. Finally, these single contact models are embedded into 
the framework of the RN in order to compute effective transport properties. All 
the proposed models—single contact as well es effective transport models—are 
verified using finite element methods. 


Regarding the pore phase, a novel method concerning the computation of the 
effective transport properties is developed. By means of the so-called Laguerre 
tessellation the pore phase of the system is decomposed into cells with each cell 
surrounding a particle. Consequently, the cell nodes and edges form the basis 
of equivalent resistor networks. The nodes are identified as pore centers and the 
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edges are the pore throats. As an extension, this model is developed further to 
account for more than one conducting species in the pore phase. Both methods 
are either verified using the finite element method or validated with the help of 
experiments taken from literature. 


It is demonstrated that the efficiency of the RN can be used to generate a large 
database with varying structure combinations. This way, the foundation is created 
for deriving prediction formulas for the effective conductivity of porous cathode 
structures represented by sphere packings with overlapping particles as well as 
the effective resistance of porous secondary particles. 


On the other hand, a mathematical model for half-cells with hierarchically 
structured cathodes with porous secondary particles is proposed. First, the 
classical half-cell model originating back to Newman and coworkers is revisited 
and the basic assumptions for the electrochemically based equations are 
presented. As a next step, the mathematical framework of the volume averaging 
method is employed to consistently extend the classical to the hierarchically 
structured half-cell model. For both models, the full set of boundary conditions 
for the half-cell setup is presented. Finally, the hierarchically structured half-cell 
model is qualitatively validated by experiments taken from literature. 


The validation of this model allows for large-scale parameter studies by varying 
electronic conductivity and diffusion coefficient of the active material as well 
as the morphology of the secondary particles. Exemplary results suggest that, 
while the rate-limiting factor for the classical cathodes is the diffusion coefficient 
of the active material, in case of the hierarchically structured cathodes, it is the 
combination of electronic conductivity and inner morphology of the secondary 
particles. 


Kurzfassung 


Hierarchisch strukturierte Aktivmaterialien in Elektroden von Lithium- 
Ionen-Batterien (LIB) mit porösen Sekundärpartikeln sind vielversprechende 
Kandidaten zur Erhöhung der gravimetrischen Energiedichte und der 
Ratenabhängigkeit der Zelle. Es gibt jedoch immer noch Aspekte dieser 
Technologie, die noch nicht vollständig verstanden sind. Um ein tieferes 
Verständnis darüber zu erlangen, wie die Kathodenstruktur und -morphologie 
die Zellperformanz beeinflusst, ist das Ziel dieser Arbeit die Entwicklung 
effizienter Werkzeuge zur Berechnung effektiver Transporteigenschaften für 
granulare Kathodenstrukturen, die in Zellmodelle importiert werden können, um 
die elektrochemische Zellleistung von LIBs zu bewerten. 


Auf der einen Seite wird die Widerstandsnetzwerkmethode (RN)—ein 
Werkzeug zur effizienten Berechnung der effektiven Transporteigenschaften 
von Partikelsystemen—hinsichtlich des Transports durch die Fest- und 
die Porenphase von granularen Medien, die durch Kugelpackungen mit 
polydisperser Größenverteilung dargestellt werden, erweitert. Was die Festphase 
anbelangt, so wird der Transport durch das Volumen der Partikel, über deren 
Oberfläche oder ein Mix aus beiden betrachtet. Für alle Fälle werden 
geeignete analytisch hergeleitete Formeln aus der Literatur verwendet oder 
entsprechend kombiniert, sodass der Widerstand zwischen zwei Einzelpartikeln 
beschrieben wird. Schließlich werden diese Einzelkontaktmodelle im Rahmen 
der RN verwendet, um effektive Transporteigenschaften zu berechnen. Alle 
vorgeschlagenen Modelle —sowohl die Einzelkontakt- als auch die effektiven 
Transportmodelle— werden mit Finite-Elemente-Methoden (FEM) verifiziert. 
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Im Hinblick auf die Porenphase wird eine neuartige Methode zur Berechnung 
der effektiven Transporteigenschaften entwickelt. Mit Hilfe der so genannten 
Laguerre-Tessellation wird die Porenphase des Systems in Zellen zerlegt, wobei 
jedes Partikel in ihnen liegt. Die Zellknoten und -kanten bilden die Grundlage 
für äquivalente Widerstandsnetzwerke. Die Knoten werden als Porenzentren 
betrachtet und die Kanten sind die Porenhälse. Als Erweiterung wird dieses 
Modell dahingehend modifiziert, sodass es möglich ist, mehr als eine leitende 
Spezies in der Porenphase zu berücksichtigen. Beide Methoden werden entweder 
mit der FEM verifiziert oder mit Hilfe von Experimenten aus der Literatur 
validiert. 


Es wird gezeigt, dass die Effizienz des RN genutzt werden kann, um eine 
große Datenbank mit unterschiedlichen Strukturkombinationen zu erzeugen. 
Auf diese Weise lassen sich erfolgreich Vorhersageformeln für den effektiven 
Widerstand von porösen Sekundärpartikeln sowie die effektive Leitfähigkeit 
von Kugelpackungen mit überlappenden Partikeln ableiten. Analog zur 
bekannten Bruggeman-Beziehung können diese Formeln in Zellmodellen 
verwendet werden, um den Einfluss der effektiven Transporteigenschaften auf 
die elektrochemische Performanz zu untersuchen. 


Auf der anderen Seite wird ein mathematisches Modell für Halbzellen mit 
hierarchisch strukturierten Kathoden vorgeschlagen. Zunächst wird das 
klassische, auf Newman und Mitarbeiter zurückgehende Halbzellenmodell 
rekapituliert und die Grundannahmen für die elektrochemisch basierten 
Gleichungen vorgestellt. In einem nächsten Schritt wird das mathematische 
Gerüst der Volumenmittelungsmethode verwendet, um das klassische auf das 
hierarchisch strukturierte Halbzellenmodell konsequent zu erweitern. Für 
beide Modelle wird der vollständige Satz von Randbedingungen für den 
Halbzellenaufbau vorgestellt. Schliesslich wird das hierarchisch strukturierte 
Halbzellenmodell durch Experimente aus der Literatur qualitativ validiert. 


Die Validierung dieses Modells ermöglicht groß angelegte Parameterstudien 
durch Variation der elektronischen Leitfähigkeit und des Diffusionskoeffizienten 
des aktiven Materials sowie der Morphologie der Sekundärpartikel. Vorläufige 
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Ergebnisse deuten darauf hin, dass, während der ratenbegrenzende Faktor bei den 
klassischen Kathoden der Diffusionskoeffizient des aktiven Materials ist, im Falle 
der hierarchisch strukturierten Kathoden es die Kombination aus elektronischer 
Leitfähigkeit und innerer Morphologie der Sekundärpartikeln ist. 
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1 Introduction 


1.1 Lithium-ion batteries 


The continuous development of numerous mobile devices like smartphones, 
laptops and tablets have a great impact on our daily life. For instance, 
almost half of the people in Germany own both a smartphone and a tablet [1]. 
We communicate with each other, schedule appointments, carry out banking 
transactions or even track our health. Overall, those little helpers have become 
an important part of our daily routine. Furthermore, fossil fuel independent 
and electrically driven automobiles continue to enter the market [2]. It can be 
expected that in some way or the other our mobility will also change. For all 
these examples lithium-ion battery (LIB) technology plays an important role. It 
is worth noting that the development and implementation of this technology has 
come a long way. 


Development of lithium-ion batteries 


Just recently, in the year 2019, the nobel prize in chemistry was awarded to John 
Goodenough, M. Stanley Whittingham and Akira Yoshino [3]. The researchers 
were honored for their contribution to the development and improvement of 
the lithium-ion battery. The committee recognized that "[t]heir research not 
only allowed for the commercial-scale manufacture of lithium-ion batteries, but 
it also has supercharged research into all sorts of new technology, including 
wind and solar power" [3]. Indeed, over the last three decades the commercial 
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and academic world experienced a large progress and a considerable amount of 
publications in that field [4]. 


Whittingham was the first who developed a battery using a titanium disulfide 
cathode (LiTiS2), a nonaqueous electrolyte and a lithium metal anode. As a 
result, the cell potential achieved was 2.5V [5]. This was an improvement 
compared to the lead-acid batteries with a cell potential of 2.0V. Later, 
Goodenough developed a cathode with an even higher potential. He used lithium 
cobalt oxide (LiCoO?) where the cell voltage increased up to 4.0V [6]. However, 
both battery compositions lacked stability due to dendrite growth of the lithium 
metal anode, which shortens battery lifetime and poses even security risks [7]. 
Finally, Yoshino used the battery setup from Goodenough, but switched the anode 
to a carbonaceous material. After running some successful safety tests, in words 
of Yoshino, this was the birth of the lithium-ion battery [8]. Indeed, modifications 
of carbonaceous materials were developed ever since [9] and the vast majority of 
anode material until today is graphite based [10]. 


Design and working principle 


There are two kinds of batteries. Primary batteries can only be discharged once 
while the secondary batteries can be recharged. Furthermore, a battery is usually 
built up of so-called cells, where they can be put in series or parallel connection 
to achieve a desired electric current or voltage. The design and working principle 
of a secondary battery cell can be observed in Figure 1.1. 


The cell has a positive and a negative porous electrode, where both of which 
consist of so-called active material. This granular material works as a host 
structure being able to take up or release lithium. Moreover, typically, 
electronically conductive carbon-black-binder mixture is added to the positive 
electrode to enhance electronic conductivity. The electrodes are separated by 
an electronically blocking porous layer, i.e. the separator. The pores of the 
electrodes as well as the separator are filled with liquid electrolyte to guarantee 
transport of the lithium ions through the whole cell. A commonly used salt in the 
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electrolyte is lithium hexafluorophosphate (LiPFg), which dissociates into Li” 
and РЕв in the solvent [10]. 


Negative Separator Positive 
electrode electrode 


Figure 1.1: Design and working principle of a lithium-ion battery cell. 


During the discharge process, lithium is oxidized at the active material surface 
of the anode. Consequently, electrons and lithium ions exit the anode material 
and move to the cathode via an external circuit and the internal electrolyte phase, 
respectively. The emerging external electric current can be used to power devices. 
Together with the oxidation process, lithium ions are reduced at the surface of 
the cathode material and the lithium atoms enter the crystal structure of the solid. 
Chemically, the above described process can be expressed as 


um tte, an 
which is, in principle, a reversible process. This means that during charging 
the above described process reverses its direction. Note that the anode and the 
cathode is defined as the electrode where the oxidation and the reduction occurs, 
respectively. Technically, only during the discharging the negative electrode 
is the anode positive electrode is the cathode. Throughout this work, only 
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discharging processes are considered which is why the anode is identical to the 
negative electrode and the cathode is the positive electrode. 


The active material, i.e. the solid particles inside the positive and negative 
electrode must be capable of taking up and releasing lithium atoms as well 
as storing them. Commonly, graphite-based materials are used in case of the 
negative electrodes [10]. Graphite is composed of multiple layers of hexagonal 
Cs structures held together by weak van der Waals forces. Lithium atoms can 
be stored between those layers where a maximum of one lithium atom can be 
stored per six carbon atoms. Therefore, the notation is „С, where 0 <x < 1 
represents the so-called degree of lithiation. 


In case of the positive electrode, the used materials are typically oxides of 
transition metals. The crystal structure must be one that allows lithium to move 
freely and be stored or released when needed. Additionally, the compound 
has to be stable for a large range of lithiation. An overview of the possible 
material structures can be found in [11]. However, new material compositions are 
continuously under development. Also, commercially used cathode materials can 
be found in [12]. For example, a commonly used material is lithium manganese 
oxide (Li; (Nii/3Mn1/3Co1/3)O» or NMC), where the capacity is high [13], the 
rate capability is good [14, 15] and it can operate at high voltages. 


A brief look into the future of batteries shows that there is further development 
in many directions. For instance, in order to tackle stability problems with the 
liquid electrolyte [16], so-called solid state electrolyte systems are a strong field 
of research [17]. Moreover, in view of the critical discussions about the resources 
used in a battery, may it be due to the abundance, geographical distribution or the 
environmental impact of lithium and other materials, the natural approach is to 
use sodium-ion batteries instead [18]. However, since those technologies are 
not market-ready yet, the chosen approach is to tweak the lithium-ion battery 
technology as far as possible to increase efficiency and thus limit resource 
consumption. 


1.1 Lithium-ion batteries 


Requirements for lithium-ion batteries 


During the past years, alternative energy sources have been pursued to reduce 
fossil fuel consumption. Therefore, electrochemical energy conversion and 
storage is under continuous development [19]. Assuming a renewable energy 
based electricity production and storage, the promising technologies are batteries, 
fuel cells and electrochemical capacitors. In this work, the focus is on lithium-ion 
batteries. 


Depending on the field of application the requirements for lithium-ion batteries 
vary. Generally, the parameters of interest are the specific energy or energy 
density in the units of Whkg~! or WhL-!, respectively. The difference between 
the two is that they are expressed either in terms of mass or in terms of 
volume. In other words, the best case scenario is having both quantities at 
high values representing high energy output while keeping weight and space 
requirements low. Moreover, these quantities can be used to compare different 
technologies against each other. For instance, the specific energy of gasoline 
reaches values of 1700Whkg~! [20] whereas current lithium-ion batteries yield 
around 250 Whkg~! and 440Whkg"! against a porous or lithium metal anode, 
respectively [20]. 


Apparently, for the lithium-ion battery technology to become a real alternative 
to, say, fossil fuels, specific energy must approximately double [21]. Since the 
specific energy is mainly governed by the positive electrode [22], finding better 
suited positive electrode materials is one solution. Those materials must have 
greater redox potentials and larger electric capacity [22]. Also, composition and 
synthesis of the material plays an important role [23, 24]. Additionally, since 
the materials remain to be poor electronic and ionic conductors, improving the 
electrode structure can lead to a better exploitation of the applied material. As an 
example, decreasing the particle size and thus increasing the surface area for the 
electrochemical reactions leads to higher rate capabilities, which means a more 
stable discharge capacity for higher currents [25]. As a drawback, however, it 
was found that cycle stability, i.e. retained capacity over multiple charging and 
discharging cycles, diminishes. Promising candidates to overcome this problem 
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are the so-called hierarchically structured cathodes where the initial material is 
processed in such a way that a distinct inner porosity of the active material is 
achieved. Primary particles, which typically built-up secondary particles of LIB 
cathodes, now form agglomerates and granular structures inside the secondary 
particles. Ultimately, in this way created cathodes experience high rate capability 
and cycle stability at the same time, see [26-30]. 


1.2 State of the art in battery modeling 


There are different methods to model batteries. On the one hand, there are 
empirical models like equivalent-circuit models. Here, simple systems of cell- 
scale ordinary-difference equations (ODEs) are used, which can be easily and 
cost-efficiently implemented in battery management systems [31]. On the other 
hand, there are physics-based models. As the name suggests, they are derived 
based on laws of physics and they allow predictions whereas the empirical 
models can only be used within experimentally marked bounds. 


Physics-based models can act on different scales [10]. The smallest of which is 
the molecular scale where the movement of atoms can be tracked and simulated. 
One scale above is the microscale or the particle-scale where the solid active 
material particles are spatially resolved, see for example [32-34]. Finally, the 
macroscale represents the cell level, where the distribution of quantities is not 
resolved spatially. Instead, the aim is to treat the cell as a continuum while 
having all inhomogeneities smeared out using appropriate volume averaging 
methods. This way, the cell can be treated as a homogeneous medium where 
the microscopic features are represented by effective transport and structure 
parameters. In this work, the macroscale modeling approach is chosen. 


1.2 State of the art in battery modeling 


Cell modeling 


A mathematical model was developed based on the works by Newman and co- 
workers [35-37], which is commonly referred to as the Newman cell model. 
Here, models for electrochemical systems were combined with porous electrodes 
using the porous electrode theory [38]. In the framework of the porous electrode 
theory, volume-averaged quantities were used where all the geometrical details 
are intrinsically accounted for [35]. Ultimately, the electrode is considered to be 
the superposition of two continua, i.e. the electrolyte and the solid phase. Later, 
the model was extended to model lithium-ion cells by [36, 39-46]. The key idea 
of the cell model is that the total electric current inside the cell is the sum of 
the ionic current, i.e. the current carried by the electrolyte, and the electronic 
current, i.e. the current carried by the solid phase. Both the currents are linked via 
electrochemical reactions at the surfaces of active material particles. The reaction 
is typically described by a charge transfer current or flux using the Butler- Volmer 
equation [47]. 


Due to the flexible nature of the cell model, it was extended in many 
directions [37]. For instance, a model to account for two different particle sizes 
within one electrode was implemented in [48] and an extension to account for 
thermal problems was presented in [49]. However, only little effort was taken to 
model hierarchically structured electrodes [50]. In [51], an impedance model for 
an agglomerate secondary particle was proposed and in [52], an electrochemical- 
mechanical model was developed for porous secondary particles. In [53], an 
agglomerate Newman type cell model was proposed. Here, an additional pore 
space was introduced on the active material particle level where the transport 
of species and the electrochemical reactions take place on the internal surfaces. 
It was possible to simulate better rate capabilities as compared to non-porous 
secondary particles. 
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Cathode structure modeling 


As mentioned above, apart from material choice, the microstructure of 
electrodes impacts cell performance [54-60]. Obviously, the particle size 
distribution influences the microstructure and thus the performance as well [61]. 
Additionally, the mechanical densification processes, such as calendering of the 
electrode, changes the microstructure [62-64]. 


Usually, by densifying the microstructure, the conductivity of the solid phase 
and therefore the electronic conductivity can be enhanced. However, this also 
reduces the pore space and thus the ionic conductivity, see [65]. In other 
words, the effective transport properties, i.e. effective conductivity and effective 
diffusivity, change with altering microstructure. Ultimately, the composition 
of the microstructure and the densification process should lead to an optimum 
balance between the effective transport properties of the solid and the pore phase. 


During the past years there have been a couple of theoretically [66-68] and 
empirically driven prediction formulas [69] to estimate the effective transport 
properties in porous media. Also, theoretical upper and lower bounds for 
effective transport properties were derived, as in [70]. Usually, the formulas take 
the form of functions of the volume fractions of the contributing phases. A well- 
known volume fraction based relation to calculate effective transport properties 
of porous media is the Bruggeman relation [71]. 


Nowadays, proper tools like the FIB/SEM tomography allow the electrode 
microstructure to be spatially resolved and virtually reconstructed. Those voxel! 
based data, the effective transport properties can be calculated using finite 
element methods (FEM) or extended finite element methods (XFEM) [72-74]. 
This way, the effective transport properties calculated are much more accurate 
than the above mentioned prediction formulas. However, the tomography, 
reconstruction and calculation process is expensive in terms of time and 


resources. 


1 Voxels are 3D extensions of 2D pixels. 


1.2 State of the art in battery modeling 


When considering the cathode structure as a granular system, a well-established 
way to model such structures is the Discrete Element Method (DEM) [75]. As 
compared to spatially resolved FEM or XFEM, in the framework of the DEM, 
each particle is represented a single element possessing a overall properties only, 
such as velocity, temperature, potential, and so on. In the most simple case, the 
elements are spheres, which are geometrically characterized by their coordinates 
and radii. Modeling granular systems as particulate assemblies is the basis of 
the Resistor Network Method (RN) [76-80]. For instance, in [76], the effective 
thermal conductivity is described by a network of contacting spheres where each 
contact is weighted by a thermal resistance. By extension, this method can 
be used for non-spherical assemblies, e.g. ellipsoids [81, 82]. As an applied 
example, in [83], resistor networks were used for the calculation of effective 
transport properties in solid oxide fuel cells (SOFC). Furthermore, in [84, 85] 
the RN was employed for the calculation of the solid phase conductivity of both 
SOFCs and LIBs. 


Resistor networks can also be used to calculate the effective transport properties 
of the pore phase [86]. Among the first was [87-89], where the permeability was 
calculated using Darcy’s law. Pore networks were created based on assemblies 
of spheres using a so-called Delaunay tessellation. The pore connecting throats 
were each weighted by a hydraulic conductivity. The hydraulic conductivity 
was determined based on the geometry based effective throat radius. A similar 
discretization approach was used in [90]. The so-called Laguerre or generalized 
Voronoi tessellation was used in [91, 92]. By means of this spatial decomposition 
technique the pore phase was divided into cells. The cell vertices and edges 
represent the pore centers and pore throats, respectively, and formed the basis of 
the resistor network. Pore throat resistances were attributed to the edges based 
on geometric bottlenecks due to the surrounding particles. 


A Statistics approach was chosen in [93, 94], where effective transport properties 
were calculated using a prediction formula based on volume fraction, tortuosity 
and constrictivity. The latter of which is a quantity resembling the influence 
of bottlenecks inside the system. To this end, numerous microstructures were 
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generated virtually and the effective transport properties were calculated by 
simulations. The generation of virtual but realistic microstructure was also 
done in [95, 96]. Also, recent developments use machine learning techniques 
on network-based modeling of the effective transport properties of particulate 
media [97, 98]. 


1.3 Objectives of this work 


The goal of this work is to combine effective transport property delivering tools 
with well-established cell models. Therefore, the resistor network method is 
extended in such a way that it becomes capable of providing all necessary 
quantities needed for cell modeling. Moreover, the Newman cell model is 
revisited and further developed in order to model hierarchically structured 
lithium-ion battery electrodes. 


The structure of this work is as follows. In Chapter 2, the fundamental 
theories and tools which are crucial for the whole subsequent work are 
presented. First, the mathematical equivalence of different transport phenomena 
is discussed. Second, the volume average method is presented, which is key to the 
development of the cell models, may it be classical or hierarchically structured 
cell models. Finally, the theoretical background of the resistor network method 
is presented and an algorithmic solving scheme is proposed. 


Chapter 3 introduces the resistor network method for both the solid and the 
pore phase. As for the solid phase, the method is extended in particular such 
that either the transport can be through the volume or via the surface of the 
particles. Additionally, acombination of both transport mechanisms is presented. 
Concerning the pore phase, the RN is applied to either the volume of the pore 
phase or to a mixture of several different species inside the pore phase. Finally, 
all the extensions are either verified using finite element methods where possible 
or validated using experimental evidence from literature. 
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1.3 Objectives of this work 


Chapter 4 is dedicated to cell modeling. The derivation of the hierarchically 
structured cell model necessitates a deeper knowledge of the classical cell model. 
Therefore, the volume average approach is used to comprehend the standard cell 
model according to Newman. Later, the same approach is used to consistently 
derive the hierarchically structured cell model. Additionally, the boundary 
conditions for both models are mathematically derived and physically motivated. 


In Chapter 5, both the resistor network method and the cell models are applied to 
specific problems. First, structural influences of secondary particles - which are 
composed of smaller primary particles - on the effective transport is studied by 
means of the RN. Second, the RN is further used to investigate the structural 
influence of densely packed sphere assemblies on effective transport via the 
particle volumes and the surfaces. Finally, the classical and hierarchically 
structured cell model is applied to real-world cathode structures. Chapter 6 
summarizes and concludes this work. 
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2 Fundamentals 


2.1 Mathematical equivalence of transport 
problems 


Under steady-state conditions, the local conservation law can generally be 
expressed by the continuity equation 


V-F=0 (2.1) 


, 


where F is the flux vector and V - (...)is the divergence operator. In case of linear 
transport, it becomes obvious that different transport phenomena can be treated 
using the identical mathematical framework, see also [68, 99]. For instance, 
energy, species and charge conservation is described by 


V-q=0, V.j—0 and Vi=0, (2.2) 


where d, j and i is the heat flux density, mass flux density and electric current 
density vector, respectively. 


The flux vectors in Equation (2.2) can be expressed via the constitutive laws 
dependant on the transport problem under consideration. The thermal, species 
and charge transport is typically defined by Fourier’s, Fick’s and Ohm’s law 


d—-AVT, j=-DVe and i=—KVo . (2.3) 
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It can be seen here that, in general, the flux is linear dependant on some 
driving force [100]. Here, A is the thermal conductivity coefficient, D is the 
diffusion, or diffusivity, coefficient and к is the electric conductivity coefficient. 
Furthermore, in the respective constitutive law, the driving force is either the 
negative gradient of temperature —VT or the negative gradient of concentration 
— Vc or the negative gradient of electric potential — V 9. 


2.2 Volume averaging method 


In this work, the focus is on porous and heterogeneous materials. In particular, 
the heterogeneities of the material are considered to be small and homogeneously 
distributed - in a statistical sense - with respect to the overall dimensions of the 
system. Itis common to refer to the smaller scale as the microscale and the larger 
scale as the macroscale. Obviously, the structure and the transport processes of 
the microscale influence the macroscopical transport properties. However, it can 
be expected that spatially resolved simulations of transport processes on the small 
heterogeneities level are computationally expensive, if possible at all. Instead, the 
idea is to subsume information of the microscale and solve the transport equations 
on the macroscale directly. Typically, a representative volume element (RVE) is 
defined where its dimensions are sufficiently large to encompass all microscopic 
phenomena and sufficiently small compared to the system geometry. Inside the 
RVEs, microscopic equations are transformed to macroscopic ones by employing 
volume averaging methods [101-105]. In the following, equations defined on the 
microscale are referred to as microscale equations whereas equations defined on 
the macroscale are referred to as macroscale equations. The structure on the 
left-hand side of Figure 2.1 represents a porous and statistically homogeneous 
material. On the right-hand side of Figure 2.1, the magnified region stands for a 
representative volume element. The macro- and the microscale is described by 
the X- and the Ё -coordinate vector, respectively, where the representative volume 
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2.2 Volume averaging method 


element is defined in terms of E . Inside the representative volume element, the 
total volume V is occupied by the &- and B-phase such that V = V^ + vB. 


a-phase 


[| 


Figure 2.1: Volume average sketch. 


Consider a general transport problem in the &-phase described by the continuity 


equation as 
9 pa 


Da +V; Fa =ba, (2.4) 


where pa = Pa (X + €,t) is the potential, Fy = Fu(&+ Et) is the flux vector, 
and ba = bg (X + E ,t) is the source term of this particular phase. The continuity 
equation must hold for every spatial point E and time f of the RVE which, in turn, 
is located at a spatial point X. 


In order to upscale the problem onto the macroscale, depicted by the x- 
coordinate, Equation (2.4) is averaged over the volume of the RVE. Practically, 
it is integrated over the volume V and divided by itself leading to 


7 / ES +v. к.) ау? = y favs, (2.5) 
or, alternatively, 
дра ay£ | ; f V fav? = 7 [ьа (2.6) 
Viv дї V Jv V Jv 
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Considering that the integral can be separated as y(...)dV® = fys(...)dV$ + 
Jyek- )dv®, where it is defined that &-phase properties are zero in every other 
phase, Equation (2.6) becomes 


1 OPa 
V Jy, Ot 


1 2 1 
dV + f| y Faavi == | Бау, (2.7) 
ү Ve V Va 


where Vg = Va (2,7) is the volume of the a-phase inside the RVE volume. 
Take note that volume averages of a property y over phase 0 are defined as 


" 1 а 
(Wa(¥,t)) = = Ya(&+&,1)dv® (2.8) 
V Ју (хи) 


апа | 
Wa = z+&,1)av®, 2.9 
Vald EN 29) 
where the former of which is called phase average and the latter of which is 
called intrinsic phase average [106]. Comparing Equations (2.8) and (2.9) it can 


be concluded that 


Val 


m 
& 
= 
| 


(Va(X,t)) = Va(X,t) = ba 5t) wa (x.t), 


(2.10) 
$a (Ху) 


where фо (X,t) is the volume fraction of the &-phase. Consequently, the phase 
averaged continuity equation in Equation (2.7) can be rewritten as 


(Bay + (V Fa) = (ba). сл!) 


As а next step, each of the terms in Equation (2.11) is transformed using volume 
averaging theorems [106, 107]. First, the first term on the left-hand side of 
Equation (2.11) is transformed into 


дра _ (Pa) 1 > 2, ё 
ae: АКГ iig dA* , (2.12) 
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where Agg = Agg(X,t) is the interfacial area between the 0- and ß-phase, 
Vag = Qr Et) and йо = rig (X 4- Et) is the velocity and the normal vector of 
the interfacial area, respectively. Note that the normal vector is pointing from 
the B- into the a-phase. The integral term represents the change of phase inside 
the RVE [106]. Assuming non-moving interfacial area and time-independent 
volume fraction of the &-phase фо, the velocity vector is equal to zero, such that 
the integral term of Equation (2.12) vanishes, yielding 


дра _ O(pa) 
Сә! д me) 


and the intrinsic volume average expression of Equation (2.13) is 


pa) _ д (Фара) a 9Pa 2.14 
ue h x 


respectively. The intrinsic volume average p, is also called the macroscopic 
potential. 


Second, the second term of the left-hand side of Equation (2.11) is converted 
using another volume average theorem, leading to 


X 1 e 
(V-Fa) =V: (Fa) — f, Fa-fadA°. (2.15) 
aß 


The integral term represents the volume-averaged flux at the interfacial area 
between the &- and B-phase. Using the average of all surface fluxes Tap» the 
integral term reduces to 


1 УТ Agp = = 
a) Fa Rada? =P Fag = aap fap. 
Aap — 


daß 


(2.16) 
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where agg is the specific surface area between the œ- and B-phase. In accordance 
to Section 2.1, the &-phase flux vector 


Fa = —kaVpa (2.17) 


is related to the gradient of the potential ро and conductivity coefficient of ky of 
phase œ. Before using another volume average theorem on the second term of 
the right-hand side of Equation (2.15), the constitutive law from Equation (2.17) 
is inserted first. Extracting the constant conductivity coefficient yields 


(кра) = —ka (Ұра). (2.18) 


The volume-averaged gradient term in the above equation is 
(Ура) = V (Pa) — "Ji Pata d4*. (2.19) 


Assuming that all particles inside the RVE have smooth surfaces and are 
symmetric with respect to their centers meaning that every normal vector on 
the surface has a counter part on the other side of the particle which points in 
the exact other direction. Also pg is assumed to be uniformly distributed across 
the surfaces. Eventually, the integral over the surface, i.e. the integral on the 
right-hand side of Equation (2.19), disappears. This leads to 


(Ура) = V (Pa) , (2.20) 


which, expressing in terms of intrinsic volume average, yields 


V(pa) = V (фаро) = aV Pa - (2.21) 


In spite of the assumptions made before, the volume fraction фо = da(X) is 
still dependent on the macroscopic location X. However, in Equation (2.21), 
the volume fraction is taken as constant and, thus, extracted from the gradient 
operator. 
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Finally, the right-hand side of Equation (2.11) simply is 
(ba) = Qaba- (2.22) 


Summarizing, the macroscopic continuity equation from Equation (2.11) is 
expressed as 


ӘР T = = 
da, +У.(—К„ф«УРе) —dopfop = Фора. (2.23) 
The term kga is defined using the effective conductivity ka eff» such that 
Equation (2.23) becomes 
Pa 2 7 b 224 
da, EN: (кара) — асв Јав = Фара · (2.24) 


The above procedure for the &-phase is performed analogously with respect to 
the B-phase, which yields 
VP pua N К И (2.25) 
08-5. + (- Ве Pp) +аав7ав = 9 В: 
In accordance to above, фв is the volume fraction, pg is ће macroscopic 


potential, kg eft 15 the effective conductivity and bg is the average source term 
of the B-phase. 


2.3 Resistor network method 


In this section, the general idea behind the resistor network method is explained 
and is partially taken from [91]. The method, as used here, follows [108], where 
the so-called node potential method is introduced. It is being used to calculate 
electric circuits. 
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From Section 2.1, irrespective of the underlying physical mechanisms, it 
becomes clear that the mathematical problems of thermal, diffusive and charge 
transport has the exact same structure. In the following, the mathematical 
formulation of the resistor network method is presented using the example 
of charge transport. When focusing for the calculation of effective transport 
properties in porous materials on charge transport, may it be electronic or ionic, 
it is understood in view of the above discussion that the method presented here 
will also apply to thermal and diffusion processes which plays an important role 
in the context of lithium-ion batteries, as well. 


Mathematical formulation of a resistor network 


Consider the example network of nodes and resistors sketched in Figure 2.2. 


9! = 90 — U9! 


4 Rx eff 


02 ^ 
D p^ 


Figure 2.2: Exemplary electrical circuit to describe the resistor network method. a) Nodes N’, № and 
the currents /^^ on the edges between those nodes. b) Potentials ф!, ф” at the nodes and 
resistances Ld at the edges. c) Using potential drop U®! = ф0 — ф! in order to calculate 
the effective transport properties. 


The nodes N’, N’ and the currents I! between those nodes, respectively, are 
indicated in Figure 2.2a. A yet unknown effective current Г. between the nodes 
№ and N! results from a potential drop between those nodes. Note that N° and 
N! have been chosen such that all the other nodes lie between them. In other 
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words, N° and N! represent the boundary nodes, i.e. the current collector nodes. 
At each node N’, Kirchhoff’s current law, 


Neigh, | 


ren Il; (2.26) 
J 


accounting for the conservation of charge, is combined with Ohm’s law, 


"E I 
Pe y N, (2.27) 


Re Кг 
representing the constitutive property of the resistors. Here, 7ineign,ı is the number 
of neighbors of node №! and U’ is the voltage between the nodes, represented 
by a potential drop ф! — @’. As a result, the equation 


Aneigh,l Aneigh,l o! = OM Nneigh,l 
"= у r= у —y-- У (v-e)6"-0 (2.28) 
7 7 Rg 7 


can be formulated for each node. Note that the directions of the currents in Figure 
2.2a can be chosen arbitrarily as long as the sign of each current in Equation 


p 


(2.28) is treated consistently. One choice would be that current pointing 


away from a node N! has a negative sign. 


According to Figure 2.2b, the unknown potentials ф', ф” correspond to the 
respective nodes М!, № and Ri! is the resistance between the nodes. Also, the 
conductance G^ = 1 IR is used as the reciprocal of the respective resistance. 
Now, it is possible to assemble a linear system of equations for the given node 
resistor network. 


In order to solve the system of linear equations, the unknown current Jer is 
replaced by an arbitrarily chosen potential drop U"! = ф0 — ọ! between the 
boundaries, see Figure 2.2c. After solving for the unknown potentials ф', the 
unknown effective current Jeg and therefore the effective resistance Rx eft = 
ud! /Tepe of the system can be calculated. 
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Solution scheme 


In the following, a general scheme for solving resistor networks according to 
above is presented. Due to its generality, the method described below is well 
suited to be used in a computer program, see [108]. 


1) Create the conductivity matrix G as 


Nneigh,l ALJ : s 
G^ if i= ; 
= : (2.29) 
—Gl otherwise 


and the current vector / as 


Tett if j=0 , 
= if yal ; (2.30) 
0 otherwise , 


where i, j — 0,1,...,n, with n being the number of nodes. Furthermore, 
Neigh, 18 the number of neighbors of the individual node. As shown in 
Figure 2.2, Jeg is the current entering the network in node 0 and leaving 
the network at node 1. The resulting system of linear equations can be set 


up as 
Goo Gio Go Спо Фо —lett 
Go. Си Gz Gni Фі let 
Go Gi? Go Gm Ф | = 0 
: : (2.31) 
Gon Gin Gon Gan Qn 0 
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2) Next, the unknown current Гов is eliminated from the right-hand side of 
Equation (2.31). As a result, a modified system of linear equations is 
created which is smaller by one degree of freedom. To this end, the equality 
ф0 = o! + U®! is used and the first row in Equation (2.31) is rewritten as 


Goo Gio Goo © Сю pı +U% —Ietf 
Gor. Gu Ga © Gm Фі Тек 
Go С; Gn > Gm Ф = 0 
Gon Gin Gan TEE Gun On 0 
(2.32) 
Second, the known voltage U°! is transferred to the right-hand side which 
leads to 
Goo Gio Goo - Gm Фі Іа — GoU’! 
Go Gu Gu cc Gm 9 Тек — Goı U! 
Go Ср Go > Gm Q |= —GoU?! 
Gon Gin Gan > Gin On —Go,U! 


(2.33) 
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Third, the unknown Kpt is eliminated from the right-hand side of the second 
line by adding the first two lines. Further, the redundant first line is deleted 


such that 
Со +Goo + Си +Gio Си Со > Gui + Gao Фі 
Соз + G12 G22 ae Gm Ф 
Gon 2 Gin Gon 22; Gan On 


e 
Se 


—(Go1 + Соо)! 
—GoU?! 


—GonU 91 
7 
(2.34) 
becomes the modified system of linear equations, where G is the modified 
conductivity matrix, @ is the modified potential vector and / is the modified 
current vector. 
3) Now, 


A 


бд=т (2.35) 
сап be solved for the modified unknown potential vector д. 


4) In the next step, the potential values have to be adjusted to the given 
boundary conditions, that is фо = 0 and фу = фо — U®!, by adding the 
offset potential Poff = Фі — фі to all ф;-о 


Фо if j=0 
g=?" | (2.36) 
Qj + Poff otherwise , 
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for j =0,...,%nodes and reassign potential values @; to Q9 where the value 
of the index j corresponds to node number J. 


5) In the last step, the unknown current к can be calculated as the sum of all 
Aneigh,i Currents entering current collector node N 0 


Nneigh 0 J 
к= У, p Б + (2.37) 
J Rx 


When the effective current Keff is known, the effective resistance Rx ef and 
therefore the effective conductance Ger can be calculated by 
y?! 1 


or Geff = 
eff Rx eff 


Rx eff = (2.38) 


Finally, as in this work the effective conductivity of a representative volume 
element is considered, domain dimensions have to be taken into account. Usually, 
rectangular domains are considered of a cross section А domain and a length Ldomain 
such that the effective conductivity can be calculated as 

Laomain legt Laomain 


Kage = С, = з 2.39 
Б i Adomain p Adomain ( ) 
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3 Modeling effective transport 
properties in assemblies of 
spheres 


A large focus of this work is the computation of transport properties in granular 
materials where the transport can be carried either by solid particles or pores. 
With modeling those materials using assemblies of spheres, the spheres represent 
the solid phase and the void space is considered to be the pore phase. In 
the following, based on the resistor network method scheme, as presented in 
Section 2.3, models are developed for the computation of effective transport 
properties in sphere packings. Since the transport can be through the solid phase 
or the pore phase, models are proposed to create equivalent resistor networks 
being applicable for both phases, accordingly. Moreover, since this is key to the 
resistor network method, appropriate models are presented to compute individual 
resistances inside the network. 


3.1 Effective transport properties 
of the solid phase 


To start with, the computation of effective transport properties of the solid phase 
is considered. Here, the transport can be through the volume or via the surface 
of the particles. Moreover, a combination of both is possible. Therefore, as 
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a first step, transport through the volume is presented. As a second step, the 
transport via the surfaces is modeled. Finally, a model is proposed representing 
a combination of both transport types. In general, the following methods use the 
geometrical bottleneck effect of contacting particles in order to calculate effective 
transport properties. 


3.1.1 Volume resistance of two 
overlapping solid spheres 


The first transport type considered is via the volume of two overlapping spheres. 
Therefore, a theory and a model of transport between two overlapping spheres is 
presented. Afterwards, this model is checked for validity using the finite element 
method [91]. 


Theory 


Many different cases of heat transfer problems were derived in [109]. In 
particular, as sketched in Figure 3.la, a steady-state heat flux фа is assumed 
through a circular aperture between two semi-infinite media. In large distance 
from the hole, a temperature Ту on the one side and 79 + AT on the other side is 
applied. In such cases, the thermal resistance can be calculated as 


E 1л + 1/42 
© 4r ; 


С. 


к, 3.1) 


The thermal bulk conductivities of the semi-infinite media are A“ and A? and 
r, is the radius of the circular aperture where the steady-state heat flux фа flows 
through. 
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3.1 Effective transport properties of the solid phase 


If, as a special case, thermal conductivity is equal to A^ = A^ = AP, 
Equation (3.1) reduces to the solution 


1 


= 
тута (3.2) 


which is given in [110]. 


Model 


As a next step, Equation (3.1) is employed to calculate the resistance of two 
overlapping spheres. Therefore, in accordance with Section 2.1, the solution of 
the thermal problem is applied to a general transport problem. In Figure 3.1b, 
the two semi-infinite media are interpreted as two half-spheres with the radii r^ 
and r^, respectively. The bulk conductivities of the /’th and J’th particle are k! 
and k”. Furthermore, the two half-spheres are geometrically overlapping and thus 
forming a contact radius r?” between them. Similar to above, a potential po is 
imposed on the middle-surface of one of the half-spheres and po 4- Ap on the 
other. The resulting flux is Fes- 


a) To 


To + AT poc Ap 


Figure 3.1: Solid-volume resistance. a) Heat flow through a circular aperture between two semi- 
infinite materials. b) Flow through two overlapping spheres. 
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Assuming that the potentials are far away from the contact radius and the volume 
of the half-spheres is comparably large, the analytical Formula in Equation (3.1) 
can be applied to two overlapping half-spheres, yielding 


ы Mg Ve 
solid,vol ^^ 4, IJ - 
Ar; 


(3.3) 
The subscript "solid, vol" refers to the transport through the whole volume of the 
contacting solid spheres. The herby described resistance shall be called solid- 
volume resistance which is based on the sphere model. 


Verification and discussion 


In order to verify Equation (3.3) and thus the solid-volume resistance, finite 
element simulations were conducted using the commercially available software 
Abaqus [111]. The FEM provides the spatially resolved solution of the stationary 
boundary value problem for two half-spheres in contact as shown in Figure 3.1b. 


The half-spheres are overlapping to form the contact radius 7. А potential 
gradient, i.e. temperature gradient AT in the FEM models, is imposed on the 
middle surfaces of the spheres. In order to calculate the resulting resistance RFFM 
in the FEM calculation, the total flux, i.e. А EM at the middle surface of one of 
the spheres is used to finally obtain 


AT 
RFEM _ 2, (3.4) 
q 


A series of finite element simulations were carried out, where not only the radius 
ratio “/-/ and the contact radius relative to the smaller sphere “/r’ was varied 
but also the bulk conductivity quotient K/W. In Figure 3.2, the dimensionless 
representation of the resistance R = т is plotted versus ће dimensionless 


contact radius “/;/” for three arbitrarily chosen cases but with relatively extreme 
ratios. 
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КРЕМ, п/п =, Kl fi = 1/1 


L1 
a 4 ve О АМ, r/r = 1/100, К/ = 1/300 
he m A REM, r/r = 1/100, K fi] = 300/1 
2 Pd —- Rig E= —0.040 


Figure 3.2: Verification of the solid-volume resistance model in Equation (3.3) with varying radii and 
conductivities using the finite element method. 


It can be seen that the values calculated according to Equation (3.3) have a 
good performance as they agree very well with the results by the finite element 
simulations. The mean error e remains within 4% for relatively large overlaps, 
e.g. 17/17 = 1, as well as for relatively small overlaps, e.g. 7/7 = 20. The latter 
of which means that the contact radius is merely 5% of the smaller of the particle 
radius. 


Furthermore, it has to be noted that the solid-volume resistance from 
Equation (3.3) even fits well with the FEM reference results when the radius 
ratio r'/r’ was varied from 1/1 to 1/100 and the conductivity ratio /&/ was varied 
from 1/1 to 1/300 and 300/1. 


In view of the above findings, it can be concluded that Equation (3.3), and thus 
the solid-volume resistance, is successfully verified to calculate the resistance 
of two overlapping spheres. It should be noted that the model makes it 
possible to not only calculate the resistance of two different-sized spheres, but 
also the bulk conductivities of the particles don't necessarily have to be equal 
anymore. Therefore, it is an improvement of the formulas proposed in, for 
example, [76, 84]. 
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3.1.2 Effective transport properties of 
overlapping solid sphere assemblies 


In the above section, it was shown that the volume resistance of two overlapping 
solid spheres can be described by Equation (3.3), i.e. the solid-volume resistance. 
For the application to an assembly of spheres, however, it is necessary that 
the formula also works well inside a network. Therefore, test cases have been 
produced and the results by RN, where the solid-volume resistance is of crucial 
importance, are compared to those by FEM analysis. Since there is no analytical 
solution of the transport problem of random spherical particle assemblies, the 
results by FEM are taken as the exact solution and thus as the reference [91]. 


Resistor network method for the solid phase 


The following method considers the solid phase of a granular structure to be an 
assembly of spheres. In contrast to, say, the finite element method, the properties 
of an individual sphere are not spatially resolved. Instead, each sphere carries 
overall properties, such as uniform temperature, electric potential, concentration 
and alike. Overlapping spheres form transport pathways through the assembly 
and, therefore, increase effective thermal and electric conductivity, diffusivity, 
etc. of the system. 


In order to calculate effective transport properties, the above mentioned transport 
- or percolated - pathways have to be converted into a network of nodes and 
resistors. As a first step, as it is highlighted in Figure 3.3a, those clusters of 
spheres have to be identified which connect the boundaries of the opposing sides. 
On the topic of cluster identification, refer, for example, to [112-114]. The grey- 
colored particles do not take part in the transport process and can be neglected. 
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Figure 3.3: Equivalent resistor network to calculate the effective conductivity of the solid volume 
phase. a) Percolating, i.e. conducting, pathways. b) Nodes, i.e. potentials, at the centers 
of the particles and edges, i.e. resistances, at the contacts, i.e. overlaps. c) Connection to 
the current collectors and boundary nodes. 


Second, as shown Figure 3.3b, the percolated clusters are converted into 
equivalent networks by assigning nodes and potentials ф', @’ to the centers of 

і А J 
the particles. Also, resistors R iiw 
nodes according to the geometric relation from Equation (3.3). Finally, as shown 


| are attributed to the edges between those 


in Figure 3.3c, additional nodes are added to model the boundary nodes where 
the boundary conditions are imposed on. The resistor network is thus created 
which is the basis for the calculation of the effective conductivity according to 
the scheme presented in Section 2.3. 


Model description 


As mentioned previously, the transport problems in this work, can be described 
using equivalent mathematical forms, see Section 2.1. Therefore, the methods 
used for, say, heat transport problems can be used equivalently for electric 
transport problems. Accordingly, the FEM solution is acquired by using the 
heat transport module in Abaqus [111], where the 10-node quadratic tetrahedron 
element DC3D10 is chosen. Similar to above, the finite element simulation 
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serves as the reference since it is the spatially resolved exact solution of the 
steady-state boundary value problem. Finally, the effective transport properties, 
achieved by both the FEM and RN methods, are compared to each other. In case 
the results by the resistor network method don’t deviate from the finite element 
solution, the former of which can be regarded as verified. 


There were two steps involved in generating the particle structures. 
First, the initial structure was created using the random close packing 
algorithm (RCP) [115]. Generally, the RCP produces randomly distributed, 
densely packed and overlap-free assemblies of spheres. Following the approach 
of [116, 117], the algorithm scope of application was extended to cover any given 
size-distribution of the particles, i.e. the radii of the spheres. In this case, a 
normal distribution with a given mean radius rmean and standard deviation ro. 


Secondly, the initial structure was further densified using an algorithm which 
will be called numerical sintering [84, 85]. While keeping the centers fixed in 
space the radii of the spheres inside the assembly were successively increased 
until a ceratin threshold was reached. The threshold can be a targeted packing 
factor, i.e. volume fraction of the solid phase ф, а, or the mean contact angle 
0, 


mean, see [85]. The contact angle 0.” of two overlapping spheres is defined 


as the bigger of the two angles 0/ and 07 enclosing the contact radius r,, see 
Figure 3.4. 


Ө” = max{07,07} 


Figure 3.4: Contact angles of two contacting spheres. 
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Finally, Ө, 


mean 18 Calculated as the mean of the contact angles of all contact pairs 


in the system. In the following, the threshold was defined using the mean contact 
angle. It is important to note that during this densification routine any sorts of 
mechanics were neglected, i.e. contact forces and alike. The geometrical data of 
the spheres were then imported into a box-shaped simulation domain and in both, 
the FEM and RN analysis, transport through two opposite conducting surfaces in 
z-direction was considered. 


Concerning the FEM models, as shown in Figure 3.6b, the spheres were cut off 
at the conducting surfaces and all nodes on one of these surfaces were set to the 
same temperature, while the potential of the nodes on the opposite surface were 
dropped by AT = 1 with respect to this potential. The resulting heat flux А EM 
was obtained at one of the surfaces where temperature was applied. By using the 
domain length Ldomain and cross section area Agomain, it is possible to calculate 
the effective thermal conductivity as 


FEM L | 
AEEM _ q domain (3.5) 


= AT Adomain 


As for the RN model, the exact same assembly of spheres can be seen in 
Figure 3.6a. In the above described framework of the resistor network method, 
the sphere assembly has to be converted into an equivalent circuit of nodes 
and resistors. In order to ensure the same electric potential at the conducting 
surfaces of the simulation domain, additional boundary nodes were created lying 
exactly on those surfaces. These boundary nodes were constructed as the contact 
points of those spheres which were overlapping with the conducting surface. 
The resistance between a boundary node and a node in the simulation domain 
1 /к! 


1,0 
Roti vot = 110 , (3.6) 


Fe 


is calculated as 


where "I? is the contact radius of the sphere J and the conducting surface, see 
Figure 3.5. 
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REO 


solid,vol 


! i! 


Figure 3.5: Construction of boundary nodes and computation of boundary solid-volume resistances. 


For the resulting network, an electric potential drop of Ag = 1 was imposed on 
the boundary nodes of two opposing surfaces. By considering the domain length 
and cross section area, the effective electric conductivity can be calculated via 


RN _ Lg Ldomain 


— “eff ботап | 3.7 
d AQ Adomain ( ) 


where the effective current Z pp was calculated using the resistor network method 
scheme from Section 2.3. 


Test cases 


For the purpose of verification, up to 15 assemblies were created, following the 
method described above. The cases were divided into certain types of assemblies, 
see Table 5.1. From the first to the third assembly type, the mean radius was 
fixed to an arbitrary length unit of 1 and the mean contact angle was fixed to 
15°, respectively. The only difference among them was the standard deviation of 
particle radii, which ranged from 0 to 0.25, such that with increasing assembly 
type number the degree polydispersity increased, as well. 
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Figure 3.6: Verification models of the solid-volume resistor network method for ап assembly of 
spheres. a) Potential distribution through the network of contacting spheres. b) 
Temperature distribution of the spatially resolved finite element mesh. 


In Figure 3.6, an example of the test cases is presented. The shown assembly 
type 3 is the one with the largest polydispersity. 


Table 3.1: Structural parameters of the spherical packings for the RN verification of the transport 
through the volume of solid spheres. 


Type Tmean Yo 0. @solid,mean Öpore,mean l'c,mean 
1 1.00 0.00 15° 0.653 0.327 0.257 
2 1.00 0.10 15° 0.650 0.350 0.245 
3 1.00 0.20 15° 0.646 0.354 0.234 
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Evaluation and discussion 


In order to compare the effective transport properties, i.e. the thermal and electric 
conductivity, provided by the FEM and the RN simulations, respectively, a so- 
called effective transport parameter К is defined. A dimensionless ratio of 
the effective transport property divided by its bulk conductivity. In case of the 
thermal transport problem and the electric transport problem, this yields 


. АЕЕМ : RN 
КРМ = апа = 0, (3.8) 
Ayulk Kpulk 


respectively. Аь is the thermal and K x is the electric bulk conductivity. In 
Figure 3.7, transport parameters resulting from FEM and RN simulations are 
plotted against each other. On the horizontal axis the results from the FEM 
solution are shown and on the vertical axis the corresponding values from the 
RN analysis are plotted. The black solid line in this figure represents a perfect 
match of both the results whereas values above or below indicate an over- or 
underestimation of the effective transport parameter by the RN with respect to 
FEM. 
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Figure 3.7: Verification of the solid-volume resistor network method using finite element simulations. 
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It can be seen that the RN results are overestimating the actual effective transport. 
However, the average error is around e = 5 % which is within an acceptable range 
considering the fact that the very complex structure is merely approximated by 
nodes and resistors. In other words, the exact solution is achieved with a lot less 
degrees of freedom and within a much shorter time [91]. 


From the above observations, a very good agreement can be stated for the resistor 
network approach with the spatially resolved finite element results. Therefore, 
the RN is considered to be verified to calculate the effective transport properties 
via the solid volume of assemblies of overlapping spheres with polydisperse size- 
distributions. 


3.1.3 Surface resistance of two 
overlapping solid spheres 


It was shown that the volume resistance of two overlapping spheres can be 
calculated by the geometrical contact radius, on the one hand, and the bulk 
conductivities of the contacting spheres, on the other hand, see Equation (3.3). 
For some applications, however, the transport doesn’t necessarily have to be via 
the volume. To increase the electronic conductivity, for example, non-conducting 
active material particles are coated with a high-conductive material [118]. In 
those cases it may be justified to assume that the transport is solely over the 
surface of the particles. 


Theory 


In [119], the electric current flow over the surface of spherical particles was 
investigated. The electric surface resistance of two overlapping spheres was 


6, 
eae Б (222 / 2) (3.9) 


2л tan(&/2) 


found to be 
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Here, Psurf is the electric surface resistivity, Ө, and 0, shall be called transport 
angle and contact angle, respectively. 


In Figure 3.8a, exemplarily, the sketch of two spheres in contact can be seen. 
It is assumed that the electric current Ý flows completely via the surface and, in 
addition, between the contact angle @, and the transport angle 6,. The transport 
angle Ө, must be in the range of [9., 180°) because otherwise the logarithmic 
expression becomes either negative or infinite. 


Model 


In order to model the transport via the surface, a suitable shell model is 
developed. The surface of a sphere is represented by a thin shell, where the bulk 
conductivity Ay, of the material under consideration is applied to. The surface 
resistivity from above is related to the bulk conductivity as 


1 
Psurf = k (3.10) 


bulk `5 
where s is the shell thickness. The idea behind this model is that as long as the 
shell is thin enough, the surface transport is accounted for. 


In Figure 3.8b, the model just described is shown. Here, two half-shells with 
the radii г! and r’ are geometrically overlapping to form the contact radius ri”. 
Now, due to the fact that the former radii can differ in size, the resulting contact 
angles 07 and 07 may also be different. The transport angles are ө; 7 and Ө; Б 
The respective thickness of the shells is given as s? and s. Also, the shell bulk 


conductivities are k’ and К”. 
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Figure 3.8: Solid-surface resistance. a) Electric current via the surface of two contacting spheres. b) 
Flow through two overlapping shells. 


Finally, the resulting resistance of two overlapping shells is calculated by a 
series connection of the surface resistance, as in Equation 3.9, where the shell- 
thicknesses and -resistances are accounted for by using the shell model from 
Equation 3.10. 


1,J I J 
Rondo = Кона surt + R solid surf 
EE tan(&”/2) ee tan(9"/2) (3.11) 
од tan(@/2) ) 2л tan(8//2) 


The subscript "solid,surf" refers to the transport via the surface of the solid 
particles. Thus, the resistance described here shall be called solid-surface 
resistance whereas the underlying model is called the shell model. 
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Verification and discussion 


A series of FEM scenarios were generated to verify Equation (3.11). Following 
the model in Figure 3.8b, the finite element models were created accordingly. 
Therefore, geometrically overlapping shells were created. To be precise, similar 
to Section 3.1.1, half-shells were created. A potential gradient was represented 
by a temperature drop of AT = 1 between the nodes at the mid-surfaces of the 
half-shells. The resulting flux is the total heat flux oF EM at one of the top or 
bottom mid-surface. It is used to calculate the resistance of the FEM models by 


AT 
КРЕМ = en: (3.12) 
q 


In Figure 3.9, the results of the investigation can be seen. Here, the dimensionless 
representation of the resistance R = Ww is plotted versus the ratio “//”, 
where г! is the smaller of the two radii. The ratio ranges from around 1 up to 
around 10, where the former means that the contact radius is almost the same as 
the smaller radius of the two and the latter means that it is 10% of the size of the 
smaller radius. 


As an example of normal to extreme cases, the radius ratio r’/r’ was varied from 
equal-sized 1/1 to relatively large ratios of !/10. Furthermore, the conductivity 
ratio ’/« was also varied from 1/1 up to 1/300. Lastly, in Figure 3.92, b and c, the 
shell thickness ratio c = s‘/r!—s! = s’/r/-s! was varied from 0.05 via 0.1 to 0.15. 
This stands for the ratio of the shell thicknesses to the inner radii which, in turn, 
are the radii of the inner surfaces of the shells. For all cases, ө” апа Ө; ” were 
set to 90°. 


LJ 
solid sur fits best for 


relatively low values of r’/r!”, that is for large overlaps. However, in Figure 3.9a, 


Generally, evaluating Figure 3.9, it can be noted that R 
it can be seen that even for a relatively large “/r/” ratio the formula performs best 


for thin shells with с = 0.05, where the mean error is around е = 3%. Also, for 
moderate shell thicknesses of с = 0.10, the mean error of e = 8 % is accaptable. 
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Figure 3.9: Verification of the solid-surface resistance model in Equation (3.11) with varying radii 
and conductivities using the finite element method. a) Shell thickness ratio c = 0.05. b) 
Shell thickness ratio с = 0.10. c) Shell thickness ratio с = 0.15. 


43 


3 Modeling effective transport properties in assemblies of spheres 


From what was found above, the solid-surface resistance and, in turn, the shell 
model is considered to be successfully verified. Using this model, it is possible to 
calculate the resistance between two thin overlapping shells, where not only the 
radii can be different, but also the conductivities of the shells. The only restriction 
is that the shells have to be thin enough. That is to say, the thickness of the shells 
should not exceed the radii of the sphere which they are covering by 15%. 


3.1.4 Effective transport properties of 
overlapping solid shell assemblies 


It was shown, that the resistance of two single contacting shells can be calculated 
using the solid-surface resistance. Next, it is investigated if the proposed formula 
in Equation 3.11 holds for an assembly of shells. 


Resistor network method for shell assemblies 


Recall that the solid phase of a granular structure was taken as collection of 
spheres. Since now the focus is on the transport via the surface of the spheres, 
rather, the structure is considered as an assembly of thin shells. Similar to above, 
overlapping shells form conducting pathways through the system and contribute 
to the effective conductivity. In order to calculate the effective conductivity of 
shell assemblies, the resistor network is used and adapted, accordingly. 


In Figure 3.10, the general approach on how to apply the RN scheme on the 
shell model is sketched. First, conducting, i.e. percolating, pathways through 
the assemblies have to be found, see the highlighted shells in Figure 3.10a. 
Next, potentials ф! and ф” are assigned to the shell centers, i.e. nodes, and the 
resistances Rs are assigned to the contact pairs, i.e. edges, as can be seen in 
Figure 3.10b. Finally, in Figure 3.10c, the boundary conditions are applied to the 
boundary nodes and the equivalent node network can be solved using the scheme 
according to Section 2.3. 
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Figure 3.10: Equivalent resistor network to calculate the effective conductivity via the solid surface. 
a) Percolating, i.e. conducting, pathways. b) Nodes, i.e. potentials, at the centers of the 
shells and edges, i.e. resistances, at the contacts, i.e. overlapping shells. c) Connection 
to the current collectors and boundary nodes. 


While using the solid-surface resistance for the contact pairs in Figure 3.10b, 
setting the transport angle between two contacting shells in the assembly needs 
special treatment. That is because the presence of more than one contact partner 
of a shell is likely to influence the transport path along the shell and decreases 
the transport angle. In Figure 3.11, exemplarily, four shells are sketched, where 
the /’th shell in the middle is in contact with the surrounding J’th, K’th and L’th 
id between the /’th and the 


solid,surf 
J'th shell, the adjacent contact pairs have to be accounted for by adjusting the 


shell. Now, in order calculate the resistance R 


transport angle ө; Jin Equation (3.11). 


Consider for the solid-shell resistance NA between shell / and J the 
contribution of the region around the J’th and Z’th shell. Apparently, the transport 
path becomes shorter, the closer the L’th shell approaches J. The length of 
the transport path is represented by the transport angle Ө, and, in particular, 


el 7L describes the transport path between the contact pair / and J and the 
neighborhood of L. Note that the transport angle ө; "TL starts from the axis 


between the centroids of /’th and J’th shell and ends at the contact radius of the 
Г and L’th shell. 
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Figure 3.11: Calculation of the transport angle for shells with adjecent contacts. 


In 3D, the angle is measured on the plane constructed by the axis between the 
centroids of /’th and J’th shell and the axis between the centroids of Pth and 
L’th shell. Since there can be more than one contacting neighbor to /, the other 
contact pairs have to be considered as well. Therefore, in a similar fashion, the 
transport angle o^ К is calculated. 


Finally, the transport angle ө; " entering Equation (3.11), is taken as the average 
of all transport angles. Thus, 
1 Nneigh 


Ө?" = — 
Ппеіеһ, 41 


0 9 for i(n) =J,K,L,...i(Meight), (8-13) 


where Mpeigh, is the number of contact neighbors of shell Z and i (n) is the index 
of the n’th shell. Note that it is chosen that Ө” " = 90° in Equation (3.13) and 


the maximim transport angle is set to 90°. 


Model description 


The generation of the shell assemblies was done in a similar fashion to 
Section 3.1.2. The RCP is used to generate randomly distributed, densely 
packed and overlap-free collections of spheres. Additionally, the particle size 
of the generated structures can be a polydisperse distribution, following a normal 
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distribution. Then, the structures were densified until a certain threshold was 
reached, in this case it was the mean contact angle. 


The densified sphere assemblies were the basis for the FEM models. In order 
to generate shells, a shell-thickness must be provided. With regard to the shell 
thickness, for every sphere imported from the generated structures, another 
sphere around the same center was imported, as well, with the shell thickness 
substracted from its radius. As a next step, the inner sphere was geometricaly 
substracted from the outer sphere and hence leaving a shell. 


In Figure 3.13b, the meshed shells can be seen. Additionally, with respect to 
the dimensions of the simulation box, the parts of the shells outside the box 
in z-direction are cut off. The color gradient in the same figure represents the 
temperature gradient. For the finite element model, a temperature gradient of 
AT = 1 was imposed on the nodes of the opposing surface of the domain in z- 
direction on the cutting plane. 


Finally, the resulting heat flux da EM. 


was evaluated at the top or bottom surface 
in Z-direction to achieve the effective thermal conductivity as 


FEM p 
FEM _ *q domain (3.14) 


eff 7 ’ 
AT Adomain 


domain length and domain cross-section area is Laomain and Adomain, respectively. 


The generated sphere assemblies from above were imported and solved using the 
RN framework. In accordance to the boundary conditions in the finite element 
model, an electric potential drop of Ap = 1 is applied to the boundary nodes in z- 
direction. The boundary nodes were generated as the contact points of the shells 
and the boundary plane in z-direction und, thus, lying exactly on that plane. The 
resistance between a boundary node and a node within the simulation domain 
was calculated as 


ro Mts tan("/2) 


К оа = 2л гап(®°/э) 


, (3.15) 
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where el? and ө; ? is the contact and transport angle between boundary plane 
and shell, respectively, see Figure 3.12. Note that ө, 0 is computed as described 
above, where the axis between the /’th shell and the boundary 0 is constructed 
between the shell center and the boundary node coordinates. 


7,0 
solid,surf 


Figure 3.12: Construction of boundary nodes and boundary solid-surface resistances. 


In Figure 3.14a, the electric potential distribution of an example assembly can be 
seen. The simulation box is depicted by the black borders. 


The resistor network analysis provided the effective electric conductivity 


IL i 

RN eff “domain 

KRN — SfE doman, (3.16) 
H Аф Adomain 


е 


where the effective current /. є was calculated using the resistor network solving 
scheme according to Section 2.3. 


Test cases 


In the following, up to 15 test cases were produced and, based on the test cases, 
the RN results were compared to the ones from FEM simulations. Three types of 
assemblies were generated. As can be seen in Table 3.3, the assemblies differ in 
their standard deviation of the radii ranging from 0.0, i.e. mono-sized particles, 
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to 0.20, i.e. polydisperse particle size distribution. In every one of the cases, the 
shell thickness ratio was set to c = s/r-s = 0.1. This corresponds to a moderate 
shell thickness according to the investigation above, which resembles the average 
case. In Figure 3.13, the RN and FEM models are compare to each other. Here, 


Table 3.3: Structural parameters of the spherical packings for the verification of the RN transport via 
the surface of solid spheres. 


Type Tmean Yo E mean Ösolid,mean @pore,mean l'c,mean Ө, mean 
1 1.00 0.00 15° 0.659 0.341 0.257 83.338° 
2 1.00 0.10 15° 0.645 0.355 0.244 84.077° 
3 1.00 0.20 15° 0.647 0.353 0.234 84.714° 


potential Z 


Figure 3.13: Verification models of the solid-surface resistor network method for an assembly of 


shells. a) Potential distribution through the network of discrete elements. b) Temperature 
distribution of the spatially resolved finite element mesh. 


49 


3 Modeling effective transport properties in assemblies of spheres 


Evaluation and discussion 


For comparison reasons, previously mentioned dimensionless effective transport 
parameters are used. Therefore, the effective thermal and electric conductivities 
are normalized as 


[FEM 2 AN d [RN KEN 3.17 
eff — 4 an eff — , (3.17) 
shell Keil 


respectively, where À pej is the thermal and куы is electric conductivity of the 
shells. 


In Figure 3.14, the results can be observerd. The effective transport parameters 
are plotted against each other. On the horizontal axis, the results from the FEM 
solution are potted whereas on the vertical axis, the RN reults are provided. The 
black solid line indicates a perfect match of both results. Points above or below 
mean over- or underestimation of the finite element results which are taken as 
reference results. 
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Figure 3.14: Verification of the solid-surface resistor network method using finite element 
simulations. 


In general, the results provided by the resistor network method overestimate the 
reference finite element results. On the other hand, the mean error ё is below 
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5%. This is a reasonably low error, given that the solid-surface model is a rather 
rough approximation of the suface transport over a very complex structure. 


3.1.5 Mixed resistance of 
overlapping core-shells 


The shell model, as presented above, can be used when the transport properties 
of the conductive coating is much larger than the core. Due to the large variety 
of different coatings for the cathode and anode material [120], for instance, 
there may be cases where the transport property of the core material can not 
be neglected [121]. 


Model 


To account for such cases, the idea now is to simply superimpose the solid- 
volume transport, Figure 3.1, and the solid-surface transport, Figure 3.8. As a 
result, in Figure 3.15, the volume part is represented by the core and the surface 
part is represented by the shell. 


Both the core and the shell radii are represented as r’ and r’. The cores have 
a conductivity of klore and k/,,. whereas the shells have a conductivity of Kon 
and ko Furthermore, the two resistances, volume and surface, are sharing the 
same contact radius r?” and contact angles Ө! and 07. The transport angles are 


ө” апа ө; 1, 
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Kren ! 


Figure 3.15: Flow through two overlapping core-shells. 


As a consequence, the superposition of the two transport phenomena leads to a 
parallel connection of the volume and the shell resistance model, Equations (3.3) 
and (3.11), respectively, yielding 


LJ E 1 1 
R lids 2 IJ + IJ 
solid,surf solid,vol 
—1 
1 | 1 
ж LJ Jl ^ VK us lo. 
khe. In ana" /2) \ , Ven In | 28/2) More Y re 
2n tan(62/2) 2n tan(6z/2) Arc 


(3.18) 
Since the equivalent resistance is а mixed form of solid-volume and solid-surface 
resistance, it will be called solid-mix resistance and the underlying model is 
called core-shell model. 
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Verification and discussion 


For verification purposes, and as was done above for the previous cases, a series 
of geometrically overlapping core-shell FEM scenarios were created. The cores 
and shells were meshed separately but shared the same nodes at the interface 
between core and shell. Hence, the temperature in the finite element simulation 
was continiuous at those interfaces. Therefore, the FEM simulation represents a 
physically more realistic situation than the superposition of the core-shell model. 


Exemplarily, in Figure 3.16, those results are shown where the radius ratios 7/7 
take moderate values, i.e. 1/1 and 1/5, and rather extreme values, i.e. 1/10. Also, 
the shell conductivities Rai and khen range both from 2 up to 10 whereas the 
cores conductivities always remain at 1. Finally, for all combinations, the shell 
thickness ratio c = 5'/4—! = s'/,/ s was varied from 0.05, via 0.1, to 0.15 to 


examine the influence of the shell thickness. 


The results of the investigation are presented in Figure 3.16. Here, the 
dimensionless representation of the resistances 


A 1 r! 
R=R 3.19 
Ven s + D m d Е УШ RE Vg as ( ) 


is plotted versus the ratio r’/r!”. R represents either the single points provided by 
the FEM simulations or the dashed lines calculated by the solid-mix resistance 
from Equation (3.18). Also, for every combination, the relative mean error e is 
shown. 


On the one hand, in Figure 3.16, it can be seen that for all cases the numerical 
results agree the best with the results by the core-shell model, if the overlap is 
the largest, i.e. r'/r!” — 1. On the other hand, even for relatively small overlaps, 
le. Mr — 10, the best agreement is achieved if the shell thickness takes its 
lowest ratio 0.05, see Figure 3.16a. For example, in case of the largest radius 
ratio “/-/ = 1/10 and the largest shell conductivities of 10, the mean relative error 
is around 146. 
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Figure 3.16: 
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as presented in Equation (3.18) 


solid,mix 


with varying radii and conductivities using the finite element method. a) Shell thickness 
ratio c — 0.05. b) Shell thickness ratio c — 0.10. c) Shell thickness ratio c — 0.15. 


3.1 Effective transport properties of the solid phase 


Influence of shell and bulk conductivity ratio 


In the following, the limits of the core-shell model, as it is used in 
Equation (3.18), are discussed. To this end, the shell to core conductivity ratio 
Keore/k pe has been varied from 1/1 up to rather extreme ratios of !/100. While on the 
other hand, the radius and shell-thickness ratios were kept constant as “/-/ = 1/1 
and c = sh fH ra 81/57 s! = 0.05. 


In Figure 3.17, the results of the above described investigation are shown. On 


top of that, the theoretical limits of the core-shell model are plotted. These 
LJ 
solid,vol 
if the volume part becomes dominant and, on the other hand, the solid-surface 
ILJ 
solid,surf 


are, on the one hand, the solid-volume resistance R from Equation (3.3), 


resistance R from Equation (3.11), if the surface part becomes dominant. 
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Figure 3.17: Influence of conductivity ratio on the solid-mix resistance model R as presented 


solid,mix 


| А А : LJ 
in Equation (3.18) and comparioson to the solid-volume model Кун vo] taken from 


Equation (3.3) as well as the solid-surface model Raus described by Equation (3.11). 


Unsurprisingly, in case of the conductivity ratio being 1, the FEM solution tends 
to the solid-volume resistance curve. Furthermore, the solid-mix resistance is 
also approaching the solid-volume curve. However, especially for small overlaps, 
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that is “//” — 10, the solid-mix resistance deviates from the FEM or the solid- 
volume resistance by around ё = 24% with respect to the numerical solution. The 
error originates from the fact that the core-shell model is a superposition of the 
sphere and the shell model. 


Another interesting observation can be made whith regard to the cases where 
the shell conductivity is 100 times larger than the core conductivity. It seems 
that this ratio is sufficient to almost exactly match the solid-surface resistance. 
Also, the relative mean error of the core-shell model to the FEM model is merely 
ё = 4%. Finally, the intermediate states are represented by the conductivity ratio 
kore/k рар = 1/10. The curve lies between the solid-volume resistance and the solid- 
surface resistance and shows a relative mean error of & = 3% with respect to the 
finite element model. 


Thus, it can be concluded that the core-shell model, i.e. the solid-mix resistance 
in Equation (3.18), is capable of covering the limits of solid-volume and solid- 
surface resistance. 


3.1.6 Effective transport properties of 
overlapping solid core-shell assemblies 


It was shown, that the resistance of two contacting core-shells can be calculated 
using Equation (3.18). Next, it is investigated if the formula holds for an 
assembly of core-shell particles. Therefore, a series of assemblies of such 
particles were generated and the effective transport properties provided by RN 
are compared to the results by FEM analysis. 


Resistor network method for core-shell assemblies 


In Figure 3.18, the resistor network approached is sketched for an assembly of 
core-shell particles. First, the percolating clusters are identified, as highlighted 
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in Figure 3.18a. The percolated clusters are then converted to equivalent node 
resistor networks. This can be seen in Figure 3.18b, where the particle centers 


аге the nodes and the contact pairs are represented by edges. Potentials ф/ and @’ 
LJ 
solid,mix 


boundary conditions are applied to the boundary nodes and the resulting resistor 


are assigned to the nodes and the resistances R to the edges. Finally, the 


network can be solved using the scheme according to Section 2.3. 
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Figure 3.18: Equivalent resistor network to calculate the effective conductivity via the solid volume 
and surface. a) Percolating, i.e. conducting, pathways. b) Nodes, i.e. potentials, at the 
centers of the core-shells and edges, i.e. resistances, at the contacts, i.e. overlapping 
core-shells. c) Connection to the current collectors and boundary nodes. 


Note that the transport angles in Equation (3.18) have to be determined in the 
same way asin the shell model from Section 3.1.4. 


Model description 


Again, similar to Section 3.1.2, randomly distributed, densely packed and 
overlap-free assemblies of spheres with a polydisperse particle-size distribution 
were generated using the RCP. Afterwards, the initial structures were further 
densified using the numerical sintering algorithm until a certain threshold was 
reached, i.e. the mean contact angle. 
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As the FEM model is concerned, the generated sphere structures were first 
imported into the simulation domain. In order to convert the spheres into core- 
shell particles, a shell thickness was provided. As a result, every sphere was 
divided into the shell and the core part. In Figure 3.20b, such a core-shell particle 
assembly can be observed. The parts of the particles outside of the borders 
in z-direction of the simulation domain were cut off. For the finite element 
calculation, a temperature gradient of AT = 1 was applied to the opposing 
surfaces of the domain in z-direction on the cutting planes. 


The resulting heat flux А EM 


direction and the effective thermal conductivity was calculated as 


was evaluated at the top or bottom surface in z- 


FEM 
_ "4 Laomain (3 20) 
Er AT Adomain і 


where the domain length and domain cross-section area is Ldomain and Адотаіп: 


The generated sphere assemblies from above were also imported to the RN 
framework. For the individual core-shell model resistance, as in Equation (3.18), 
the shell thickness has to be provided. Similar to the FEM model, the potential 
drop Ap = 1 was applied to the boundary nodes of the opposing surfaces in z- 
direction. The boundary nodes were generated as the contact points of the core- 
shells particles and the boundary plane in z-direction und, thus, lying exactly 
on that plane. The resistance between a boundary node and a node within the 
simulation domain was calculated as a parallel connection, again, according to 


—1 


1 
1,0 
Ro = | (3.21) 
RER Val s ln tan(6/ 9/2) | Vikore 
2л tan(02/2) CA 10 


where o1? and ө; ? is the contact and transport angle between boundary plane 
and shell, respectively, and 19 is the contact radius, see Figure 3.19. 
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1,0 
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Figure 3.19: Construction of boundary node and boundary resistance. 


Finally, the RN provides the effective electric conductivity using 


I Ldomain 
Kore = lt (3.22) 
TET AQ A domain 


where ^+ was calculated using the resistor network solving scheme from 
Section 2.3. 


Test cases 


Up to 15 test cases were created and, based on the test cases, the RN results were 
compared to the ones from FEM simulations. Three types of assemblies were 
generated. As can be seen in Table 3.5, the assemblies differ in the standard 
deviation of the radii ranging from 0.0, i.e. mono-sized particles, to 0.20, i.e. 
polydisperse particle sizes. In every one of the cases, the shell thickness ratio 
was taken as с = s/r-s = 0.1 representing a moderate shell thickness from the 
above investigation. In Figure 3.20, the RN model is compared to the same by 
the FEM model. Here, the assembly is of the type 3 with the largest standard 
deviation rg. 
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Table 3.5: Structural parameters of the spherical packings for the verification of the RN transport 
through the solid core-shells. 


Type "mean ro B mean Ösolid,mean Öpore,mean lc. mean Ө, mean 
1 1.00 0.00 15° 0.655 0.344 0.257 83.35° 
2 1.00 0.10 15° 0.647 0.353 0.242 84.10° 
3 1.00 0.20 15° 0.647 0.353 0.233 84.74° 


potential Z 


Figure 3.20: Verification models of the solid-mix resistor network method for an assembly of core- 
shell particles. a) Potential distribution through the network of an assembly of core-shell 
particles. b) Temperature distribution of the spatially resolved finite element mesh. 
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Evaluation and discussion 


Before comparing the effective transport properties, they are normalized by th 
bulk conductivities first. In this case, the bulk conductivities were taken from the 
core parts. This yields the dimensionless effective transport parameters as 


;FEM Adr “RN KN 
_ Te _ Хе 
Keft,mix m À and keff mix — E (3.23) 
core core 


respectively, where A ore is the thermal and Kore is electric conductivity of the 
cores. The results of the above investigation can be seen in Figure 3.21. The 
effective conductivity is overestimated by the RN. However, the mean error is 
around e = 7% which is acceptable given the fact that the core-shell model is 
a rough simplification. In this sence, the core-shell model can be regarded as 


verified to a degree sufficient of this work. 


0.875 0.900 0.925 0.950 0.975 1.000 1.025 1.050 


[FEM 
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Figure 3.21: Verification of the solid-mix resistor network method using finite element simulations. 
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3.2 Effective transport properties 
ofthe pore phase 


In this section, the effective transport properties of the pore phase are 
investigated. The pore phase is considered to be the space between the solid 
phase, e.g. active material particles in lithium-ion battery cathodes. Here, 
this space is filled with liquid electrolyte which is needed to provide sufficient 
lithium-ions and ionic charge for the electrochemical reactions on the active 
material particle surfaces. In other words, the electrolyte-filled pore phase has 
to provide a sufficiently large effective diffusivity and conductivity for a lithium- 
ion battery to operate properly. 


In addition to the electrolyte, inside the pore phase an electronically conductive 
carbon-black-binder phase, i.e. filler phase, is distributed. This is in particular 
necessary since the active material particles are commonly known for their poor 
electronic conductivity. This, on the other hand, is needed to provide enough 
electrons for the electrochemical reactions. The pore phase is thus shared by 
both the liquid electrolyte and the filler phase. 


In the following, a resistor network method is developed for the calculation of the 
effective conductivity of the pore phase. Additionally, an extension is proposed 
to account for those cases where the pore phase is shared by multiple conducting 
phases. 


3.2.1 Volume resistance of pore throats 


The crucial part of the resistor network approach is the identification and 
computation of individual resistances. In the case of the pore phase, the 
individual resistances are the pore throats. Therefore, a model is proposed to 
estimate the resistance of such pore throats based on the geometrical bottleneck 
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effect [91, 92]. For simplicity, the model is presented for the 2D case. The steps 
necessary for a 3D implementation are mentioned whenever needed. 


Model 


To start from, inside the pore phase of a granular assembly, it is assumed that 
transport is hindered by the bottlenecks created by the solid particles. To model 
such bottlenecks, i.e. pore throats, it is necessary to identify the pores and pore 
centers first. 


a) b) с) 
F 
N! N’ 
Fres 
AL” Amean 
Pore 

throat Sub- 

throat 

N! J N! 
R reval 


Figure 3.22: Pore-volume resistance. a) Flux through a bottleneck between two particles. b) 
Resulting flux between nodes №! and N” inside the spatially discretized bottleneck, i.e. 
pore throat. c) Discretization of throat into sub-throats. Divide sub-throats into small 
increments in order to calculate wire resistances. 


Figure 3.22a shows an example of such a bottleneck where two particles form 
a narrow region, i.e. a pore throat, between them. Before and after the pore 
throat are the pore centers which are connected by the pore throat. Between the 
pore centers, some flux F is considered to flow through. The exact position of 
the pore centers is achieved by a spatial decomposition technique, which will be 
subject to the next part. Next, in Figure 3.22b, the nodes N’ and N” indicate the 
start and end point of the flux. Those nodes are taken as the pore centers and the 
gray-shaded region is the spatial decomposition of the pore throat. 
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A pore throat resistance Bay is attributed to the pore throat which is 
expressed as a parallel connection of resistances representing the sides of the 
node connecting edge. In the 2D example, see Figure 3.22b, these sides are the 
sub-regions created by triangles of the edge defining nodes and the respective 
centers of the surrounding particles. Additionally, the overlapping areas of the 
particles are subtracted from the triangles. Because the sub-regions describe parts 
of the pore throats, they are called sub-throats. 


A sub-throat is sketched in Figure 3.22c. The resistance of such a sub-throat is 
calculated as a series connection of wire resistances. To this end, each edge is 
divided into sufficiently small increments AL” and the resistance АК” of each 
increment is calculated by 
AL” 
АК" = Ppore Zn — ; (3.24) 


mean 


where Ppore is the resistivity of the pore material. A mean line length Ат of 
the n’th increment has to be passed, in case of 2D. In case of 3D, its a mean cross 


section area. 


А А LJ i s; 
The resulting resistance Ау; of one sub-throat m is then calculated as a series 
connection of the incremental resistances as 


RIJ = Nincr AR” Е Niner AL” His 
„=, m = Ppore У, An , (3.25) 
n=1 = 


where nincr is the number of increments between nodes N’ and N”. Here, AL is 
estimated individually for each sub-throat. The size of AL is successively reduced 
until the resulting resistance calculated for the next increment doesn't differ from 
the previous one by a certain small threshold, here within 5 %. 


Finally, as mentioned above, the resulting resistance 


-1 


TJ Msthr 1 
Кооте,уо1 = Y EJ (3.26) 
m1 Rm 
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of the edge, i.e. pore throat, which connects the nodes М!, №, ie. the pore 
centers, is a parallel connection of the resistances Ri of all тим sub-throats 
belonging to this edge. Thus, a model is derived for pore throat resistances which 
can be used along the resistor network method when modeling the pore phase 
effective transport properties. The index "pore,vol" refers to transport through 
the volume of the pore phase which is why this resistance is called pore-volume 
resistance. 


3.2.2 Effective transport properties of 
pore throat assemblies 


In the above section, the computation of an individual throat resistance was 
presented, i.e. the pore-volume resistance from Equation (3.26). Next, this kind 
of resistance is used in the framework of the resistor network method to calculate 
effective transport properties of the pore phase. 


Resistor network method for the pore phase 


While it is quite straight forward to discretize the solid phase of a granular system 
using spheres, as demonstrated in Section 3.1, it is not as obvious to model the 
pore phase in a similar way, because of the relatively complex shapes of the pores 
and throats. In Figure 3.23, exemplarily, an equivalent resistor network for the 
pore phase is sketched. 


First, pores and pore connecting throats have to be identified. To this end, the 
pore phase is discretized using the so-called Laguerre or generalized Voronoi 
tessellation method [122, 123], where the generators are given by the midpoints 
of the particles and their weights are given by the corresponding radii. For 
detailed information regarding this tessellation technique, refer to [124]. The 
computation is done by employing the software library Voro++ [125]. 
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Figure 3.23: Equivalent resistor network to calculate the effective conductivity of the pore volume 
phase. a) Discretization of the pore phase using Laguerre cells. b) Potentials at 
cell vertices and resistances at the edges. c) Connection to the current collectors and 
boundary nodes. 


In general, for isolated points in a domain, which happen in this case to be 
the sphere centers, the tessellation method assigns a spatial cell to each such 
point. Each cell contains all points whose distance to the associated cell center 
is less than or equal to any other sphere or cell center. In case of packings with 
polydisperse particle sizes, additionally, the distances are adjusted with respect to 
the particle sizes [122]. In 3D space, the cells take the form of convex polyhedra. 
Using this method, the whole domain volume can be fully decomposed into cells. 
Finally, the pore phase is defined as the volume of the cells where the volume of 
the particles, i.e. solid phase, is subtracted from. 


The result of the tessellation is that all particles are wrapped in cells while their 
borders, i.e. edges and faces, lie in an optimum distance between the cell centers, 
i.e. particle centers. An exemplary cell discretization can be seen in Figure 3.23a. 
The spherical particles are wrapped by cells, edges, vertices and, in the 3D-case, 
faces. Each vertex is chosen as the center of a pore and the corresponding edges 
as the pore throats. By extension, the network nodes N! and № correspond to the 
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vertices where the potentials g/, ф” are associated with and the edges are related 


А LJ А 
to the resistances Rare var see Figure 3.23b. 


In order to compute the resistances using the approach described above, the 
regions around the individual edges have to be decomposed into throats and sub- 
throats. In case of 2D, the regions are constructed as the areas which are given by 
triangles defined by the nodes of the considered edge and the associated sphere 
centers, where the intersection area with the associated spheres is subtracted 
from. In the 3D case, the surrounding region is the volume given by tetrahedrons, 
which, as an extension of the 2D case, are additionally defined by the centers of 
the faces of the cells which meet at the considered edge. The intersecting volumes 
of the associated spheres are then subtracted from the tetrahedrons. 


The boundary nodes have to be identified where the boundary conditions can be 
applied to, see Figure 3.23c. Finally, the resulting node resistor network is solved 
by the scheme described in Section 2.3. 


Verification 


As a next step, the verification of the pore-volume model in the framework of 
the resistor network method is checked. To this end, similar to Section 3.1.2, 
different assemblies of varying degree of particle polydispersity were generated 
using the RCP and densified using the numerical sintering algorithm. The 
resulting collections of spheres were then imported to both the RN and the FEM 
framework. Finally, the resulting effective transport properties were computed 
and compared to each other. 


As far as the finite element model is concerned, first, the generated sphere 
assemblies were imported into a box-shaped domain. Secondly, the spheres were 
cut out of the domain such that the complement volume remained. The latter of 
which was identified as the pore phase and can be seen in Figure 3.24b. The color 
gradient refers to the temperature distribution emerging from the temperature 
drop of AT = 1 on the surface nodes of the simulation domain in the z-direction 
while the other surfaces in x- and y-direction remain insulated. The resulting 
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heat flux ól EM 


was evaluated at the top or bottom surface in z-direction and the 
effective thermal conductivity was calculated as 


FEM р 
АЕЕМ _ 7q domain (3.27) 


eff T , 
AT A domain 


where the domain length and domain cross section area is Laomain and Aaomain; 
respectively. 


Regarding the resistor network model, the spheres were imported into the box- 
shaped simulation domain, as well. As a next step, the Laguerre tessellation was 
used, as described above, to fully discretize the pore phase and which was the 
basis for the resistor network method. In Figure 3.24a, the potential distribution 
of such a resistor network can be observed. The vertices on the boundary planes 
in Z-direction were directly used as boundary nodes and an electric potential drop 
of Ag = 1 was applied to these boundary nodes. 


The effective electric conductivity was calculated by 


FEM . La Laomain 


K, = — —_——_, 3.28 
= АФ Adomain ( ) 


where the effective current Z pp was calculated using the resistor network solving 
scheme according to Section 2.3. 


Test cases 
Next, 15 test cases were created, where the structural parameters upon which 


the assemblies were generated are presented in Table 3.7. In Figure 3.24, as an 
example, one assembly of the type 3 with the largest standard deviation is shown. 
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Table 3.7: Structural parameters of the spherical packings for the verification of the transport through 
the volume of the pore phase. 


Type "mean Yo Omean @solid,mean Öpore,mean 
1 1.00 0.00 15° 0.657 0.343 
1.00 0.10 15° 0.649 0.351 
3 1.00 0.20 15° 0.644 0.356 


potential Z 


Figure 3.24: Verification models of the pore-volume resistor network method for the pore phase of an 
assembly of spheres. a) Potential distribution through the network of pores and throats. 
b) Temperature distribution of the spatially resolved finite element mesh. 


Evaluation and discussion 


In order to compare the RN and FEM results, the effective transport properties 
have to be converted to effective transport parameters. Therefore, the effective 
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thermal and electric conductivity were normalized by the bulk thermal А, ц, and 
electric conductivity къу respectively, which yields 


А А. FEM ^ КРЕМ 
КРМ = апа =. (3.29) 
ulk Kbulk 


In Figure 3.25, the comparison of the investigation can be seen. As usual, the 
finite element results are on the horizontal whereas the resistor network values 
are on the vertical axis. The solid black line indicates a perfect match of the two 
methods. Here, the results obtained by the RN are overestimating the reference 
FEM results by a little bit as the mean error is around е = 1%. This, on the 
other hand, is a low error considering the approximation of the arguably complex 
pore phase structure by a combination of series and parallel connections of wire 
resistances. Moreover, concerning resource and time cost, the RN is superior to 
the FEM calculation [91]. 


0.220 0222 0224 0.226 0.228 0.230 0.232 0.234 0.236 
{РЕМ 


Figure 3.25: Verification of the pore-volume resistor network method using finite element 
simulations. 


Given the rather low error, even for the pore phase of polydispersely sized sphere 
packings, the resistor network method to calculate the effective conductivity of 
the pore phase can be considered as verified. 
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3.2.3 Volume resistance of mixed pore throat 


Different other material phases can distributed inside the pore phase. For 
example, to enhance the electronic conductivity and to fix the active materials 
of the cathodes of lithium-ion batteries, solid carbon-black (CB) and polymeric 
binder is added to the pore phase. Thus, the pore phase is shared by liquid 
electrolyte, on the one hand, and by a conductive filler phase, i.e. carbon-black- 
binder phase, on the other hand. In the following, the model of the pore-volume 
resistance is extended in such a way that it is able to account for such cases. It is 
assumed that the different phases in the pore phase don’t interact with each other. 


Model 


In Figure 3.26, the bottlenecks, similar to Section 3.2.1, are sketched. The only 
difference is that another phase is added to the bottleneck space. The added 
phase is indicated by the black dots in Figure 3.26a. Next, in Figure 3.26b, the 
pore throat is identified as the gray-shaded area. 


a) b) c) 
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Figure 3.26: Pore-mix resistance. a) Flux through a bottleneck between two particles. Additional 
phase depicted as black dots. b) Resulting flux between nodes N’ and № inside the 
spatially discretized bottleneck, i.e. pore throat. c) Discretization of throat into sub- 
throats. Divide sub-throats into small increments in order to calculate wire resistances. 
Account for volume fraction of the additional phase. 
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LJ 
pore,mix,c? 
refers to the transport through the particular phase c of a mixed pore phase. In 


The resistance of such a pore throat with added phases is R where c 
order to calculate the resistance, again, each edge is divided into sufficiently 
small increments AL” where an incremental resistance AR; of each segment 
is calculated by 


AL” 
N,c с 1 
АК = Ppore® Zn — + (3.30) 
mean 
Here, Ppore is the resistivity of the phase c, 67, is the volume fraction of the phase 
inside the sub-throat m and Ал, is the mean cross section area for the n'th 
increment. 


The resulting resistance Rise of the added phase c inside the sub-throat m is then 


calculated as a series connection of the incremental resistances as 


1Ј Nincr Nincr n 

at ME AC e 

Riv = У АЕ = Ppore? У, а, (3.31) 
n=1 n—]| ^^mean 


where the number of increments used between Nodes №! and №” is niner and AL 
is chosen individually for each sub-throat. For this, AL is successively reduced 
until the resulting resistance of the sub-throat doesn't differ from the previous 
one by a certain small threshold. In this work, the threshold is set to 5 96. 


LJ 


The resulting resistance R reine 


of the added phase c inside the pore throat is a 


parallel connection of the resistances RI/« of all Msthr sub-throats. 


—1 


TJ "г 1 
R pore.mix,c = у 1,7, 5 (3.32) 
т=1 Ёт 


Summarizing, for each phase of the роге phase separately, а роге throat resistance 
was derived where the transport path is shared by more than one phases. The 
index "pore,mix" refers to transport through the mixed volume of the pore phase 
which is why this resistance is called pore-mix resistance. 
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3.2.4 Effective transport properties of 
mixed pore throat assemblies 


In the following, the pore-mix resistance is used inside the framework of the 
RN. To this end, a pore-mixed resistor network is presented and validated using 
experimental results from the literature. 


Resistor network for the pore mix resistance 


The build up of the resistor network for the mixed pore structure is similar to 
Section 3.2.2. As can be seen in Figure 3.27a, the Laguerre tessellation is 
used to discretize the pore phase. Recall that the resulting cells are formed by 
the vertices, edges and faces. The black dots represent an additional randomly 
distributed phase inside the pore phase. 
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Figure 3.27: Equivalent resistor network to calculate the effective conductivity of the pore mix phase. 
a) Discretization of the pore phase using Laguerre cells. b) Potentials at the Laguerre 
vertices and representative resistances at the edges. c) Account for volume fraction of 
added phase in each sub-throats. Connection to the current collectors and boundary 
nodes. 
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In Figure 3.27b, it can be seen that the vertices are taken as the nodes where the 
potentials ф', Ф” are associated with. The edges stand for the representative pore- 
eee Its computation is performed using the methodology 
described in Section 3.2.3. To this end, the distributed additional phase inside 
the pores is regarded as smeared across the sub-throats. This is depicted by 
the differently shaded areas, i.e. sub-throats, in Figure 3.27c. Note that the 


differently shaded areas imply that, if known, each sub-throat can have individual 


mix resistance R 


volume fractions of the considered phases. Finally, the boundary conditions can 
be applied to the boundary nodes and the effective conductivity can be calculated 
using the scheme described in Section 2.3. 


Validation 


In the above sections, the FEM was rigorously used to verify the individual 
resistor network approach. However, in case of lithium-ion batteries, this 
is not feasible due to the very complex structure of the mixed phases in 
the pore phase. Alternatively, the resistor network approach for the mixed 
pore resistance is validated using experiments from literature. In [126], 
the effective electric conductivity of carbon-black-binder mixtures with and 
without the presence of active material particles was measured. Several binder 
compositions were presented. In this case, the binder consisted of polyethylene 
oxide (PEO), polyvinylidene difluoride-co-hexafluoroproylene (PVdF-HFP) and 
ethylene carbonate-propylene carbonate (EC-PC). The carbon-black content of 
the mixtures was varied from 0 to 18.8%. 


It can be seen that for a carbon-black volume fraction of 0%, the effective electric 
conductivity does not vanish. That is because the ionic conductivity of the binder 


-l. Since 


contributes to the electric conductivity, which is around 1-10~7Sm 
this compositions have no relevance in applications, in this investigation, carbon- 
black-binder mixtures and electrode compositions without carbon-black content 
are neglected. Therefore, if the carbon-black volume fraction reaches zero, the 
electric conductivity is also set as zero. This way, the electronic conductivity is 


extracted from the experimental data. The effective electric conductivity rises 
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orders of magnitudes for carbon-black volume fractions from 1.5 to 6%. This 
range was identified as the so-called percolation threshold, which resembles 
the critical volume fraction of conducting species where the likelihood to form 
conducting pathways increases drastically. For contents larger than that, the 
effective electric conductivity increases asymptotically until, at the theoretical 
limit of carbon-black volume fraction of 100%, the bulk electronic conductivity 
of carbon-black is reached, which is assumed as 1000S т! [126]. 


In order to account for the above described behavior, the effective electronic 
conductivity of the carbon-black-binder mixture provided by the experimental 
data was fitted by the s-shaped function 


KCp-B.ett =a: | 1— exp (- (=) | , (3.33) 


where the parameters аге a = 1000, b = 0.167 and с = 2.358. Here, a was chosen 
as the bulk electronic conductivity of carbon-black. Thus, for large values of 
the carbon-black volume fraction фсв related to the carbon-black-binder mixture 
the bulk electronic conductivity would be reached. The fit-function fulfills the 
criterion that the effective electric conductivity is OS m^ if the carbon-black 
content reaches zero, see Figure 3.28b. 


As the next step, active material particles were added to the carbon-black-binder 
mixtures and the effective electric conductivity was measured. In the following, 
the electrode compositions are divided into two phases, which are the pore and 
the solid phase. The former of which comprises of the electrolyte and the carbon- 
black-binder mixture phase and the latter of which is the active material phase. 
It was reported that the electrode samples were densified such that porosities of 
30% were achieved. Since the pores are filled with electrolyte, the resulting 
phase is identified as the electrolyte phase. The remaining 70% of volume 
fraction were shared by the carbon-black-binder mixture and the active material 
phase. This phases contained a volume fraction of 19 — 22% active material 
and 81 — 78% of carbon-black-binder mixture. For the following investigation, 
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the active material and carbon-black-binder volume fraction with respect to the 
whole volume was set as the average values, which were around 15% and 55%, 
respectively. 


Experimentally achieved effective electric conductivity values of the above 
described electrode structures were compared to numerical results. To this 
end, comparable virtual compositions were created. Similar to Section 3.1.2, 
the initial sphere packings were generated using the RCP. The structures were 
further densified by means of the numerical sintering algorithm. Additionally, 
the algorithm was slightly modified. In this case, the box-shaped simulation 
domain was isotropically and successively reduced or expanded such that the 
spheres inside the domain were moved upon or apart from each other. This was 
done until the targeted porosity was achieved. This way, the particle size was 
preserved. 


The used active material was Li; 2 V3Og and the size was set to a particle radius of 
lum, which was deduced from SEM images in [126]. Following the description 
of the electrode composition from above, the targeted volume fraction of the 
active material with respect to the whole volume was set to 15%. The remaining 
pore phase of 85% was shared by 65% and 35% carbon-black-binder mixture 
and electrolyte, respectively. 


Based on this composition, scenarios of sphere assemblies were created which 
were randomly distributed but structure-wise identical to the electrode in the 
experiments. Using the pore-mix resistor network method, the presence of 
added phases, i.e. electrolyte and carbon-black-mixture phase, was accounted 
for by appropriately setting the volume fraction c in the pore throat resistance 
computation from Equation (3.26). In order to calculate the (effective) 
bulk conductivity of the carbon-black-binder phase, the fit-function from 
Equation (3.33) was used. 
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Figure 3.28: Validation of the pore-mix resistor network method by using experiments from [126]. 
a) Effective electric conductivity of a cathode structure with active material and carbon- 
black-binder phase. b) Effective electric conductivity of the carbon-black-binder phase 
as a function of the carbon-black content, see Equation (3.33). 


In Figure 3.28a, the results of the resistor network calculations, i.e. KS св-в, are 
compared to the measurements from [126]. It can be seen that the values provided 
by RN match the experiments very good for carbon-black volume fractions 
above 4%, which is within the observed percolation threshold [126]. However, 
the numerical results deviate a little bit for carbon-black contents below 4%. 
This originates from the fact that the fit-function used for the conductivity of 
the carbon-black-binder phase deviates for lower carbon-black volume fractions. 
On the one hand, the fit-function was chosen such that it represents the s-shaped 
percolation phenomenon, where it is difficult to properly describe the behavior at 
the vicinity of the percolation threshold. On the other hand, usually, carbon-black 
volume fraction are above 4%, see for instance [127]. In general, the overall 
quality of the match between the numerical results and the experiments leads 
to the conclusion that the resistor network approach described above is able to 
calculate the effective conductivity of added conductive phases in the pore phase. 
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3.3 Summary 


In Sections 3.1 and 3.2, it was shown how the resistor network method can be 
used to model effective transport properties of both the solid and pore phase in 
an assembly of spheres with polydisperse size-distribution. 


Concerning the solid phase, equivalent resistor networks were created based on 
boundary connecting clusters of spheres. The transport through this clusters 
was distinguished between volume and surface type. In order to be used in 
the framework of the RN, the crucial part is the computation of resistances of 
individual contact pairs. To this end, regarding the volume transport, the solid- 
volume resistance was constructed in Equation (3.3) based on the sphere model. 
The surface transport was computed using the solid-surface resistance from 
Equation (3.11) where the underlying model was the shell model. Moreover, 
assuming the transport is partially through the volume and the surface, a solid- 
mix resistance was developed in Equation (3.18). The corresponding model was 
the core-shell model. 


As for the pore phase, equivalent resistor networks were created based on spatial 
decomposition techniques. Here, the Laguerre tessellation was used where every 
particle inside an assembly is placed into convex polyhedral cells. This way, 
the pore phase was geometrically defined as the complement volume of cells 
and the contained spheres. Additionally, the cell vertices and edges were the 
basis upon which resistor networks were created. While the vertices resembled 
the pore centers, the edges were identified as the pore throats. Based on the pore 
throat geometry, aresistance was calculated, accordingly. Since the transport was 
through the volume of the pore, the corresponding resistance in Equation (3.26) 
was named pore-volume resistance. Additionally, since the pore phase can be 
shared by multiple conducting phases, another model was proposed. The pore- 
mix resistance in Equation (3.32) considers the resistance of additional phases 
using their volume fractions and resistivity. 
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The working principle of lithium-ion battery cells was described in Section 1.1. 
In general, electrochemical reactions take place on the surfaces of the active 
material of the positive and negative electrode. Those electrodes usually are 
porous structures and comprise of different components. In order to model 
battery cells, several physics-based approaches on different length-scales can be 
applied [10]. In the following, a continuum approach is chosen. 


A prominent candidate of this modeling approach is the Newman type cell 
model [41]. While the development of novel cathode structures proceeds, i.e. 
hierarchically structured electrodes [26-29], the need for novel cell modeling 
approaches arises, as well. Therefore, in the first section of the following part, the 
classical cell model is presented first. Next, In the second section, the cell model 
is consistently extended to account for hierarchically structured electrodes. 


4.1 Classical half-cell 


In Figure 4.1a, the sketch of the classical half-cell setup is shown. The anode 
on the left-hand side is a lithium foil and the positive electrode, i.e. cathode 
during the discharging process, on the right-hand side is a porous composite 
structure, comprising the active material and carbon-black-binder mixture, which 
is the conductive filler material. The porous separator is sandwiched between 
the lithium foil and the positive electrode. In addition, both the porous cathode 
and separator are completely soaked with liquid electrolyte. The setup is chosen 
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such that one is able to study solely the influence of the porous cathode on the 
performance of the cell. In contrast, the full-cell setup, with an additional porous 
anode, both electrodes would influence the cell performance. 


WIDER 


Anode Separator Cathode Cell 
level 


Particle 
level 


Figure 4.1: Half-cell setup of lithium-ion batteries. a) Sandwich structure of a lithium-ion battery cell. 
The anode is a solid lithium metal. b) The porous structure of the cathode is homogenized. 


A suitable model to incorporate battery charging and discharging processes 
originated from [36, 40-42]. The model basis is the so-called porous electrode 
theory [47] and [35]. In this theory, the structural details of the actual geometry 
is smeared out across the model to achieve effective properties. This way, two 
scales are introduced, where transport defined on the full three dimensional 
structural details of the electrode on the smaller scale is called microscopic while 
the corresponding effective transport on the larger smeared scale is referred to 
as macroscopic. Therefore, transport equations defined on the microscale are 
referred to as microscale equations whereas transport equations defined on the 
macroscale are referred to as macroscale equations. Every point of the model 
represents a superposition of two phases. While the electrolyte represents one 
phase by itself, typically, the active material and filler material is subsumed as 
the other phase, i.e. the solid phase. The structural and material properties of 
each phase are characterized by the volume fraction, specific surface area and 
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effective transport properties. Additionally, the two phases interact with each 
other via electrochemical reactions on the active material surfaces. 


Thus, the particle and pore geometry of the system is not spatially resolved 
but, rather, can be viewed as smeared across the model. In Figure 4.2b, the 
smeared region is indicated by the shaded area where, additionally, at every 
point in the cathode region the active material particles are modeled as spheres. 
Therefore, the cell model, as presented here, is divided into two levels, i.e. the 
cell and the particle level. The former of which represents the macroscopic scale, 
i.e. transport in the porous electrode, whereas the latter describes the transport 
processes inside the active material. 


4.1.1 Macroscale equations 


In the framework of the porous electrode theory, at every spatial point of the 
macroscale model, the solid and the electrolyte phase are superimposed. By 
extension, the currents and fluxes in the phases are superimposed, as well. It is 
assumed that, on cell level, the electronic current is carried by the solid phase, 
whereas the ionic current is carried by the electrolyte phase. Furthermore, while 
the transport of cations is supposedly via the electrolyte phase on cell level, the 
insertion process and lithium transport process inside the solid phase is modeled 
on particle level. 


The above transport phenomena are described by four continuity equations 
representing the conservation of electronic and ionic charge as well as the 
conservation of mass in the solid and electrolyte phase. In the following, the 
micro- and macroscale forms of the continuity equations are recalled as originally 
presented in [35, 41, 47]. As it was shown in [128, 129], the macroscopic forms 
used in the cell models can be derived by applying volume averaging methods, 
which were presented in Section 2.2. The same methods are also chosen in the 
following. 
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Conservation of mass in the solid phase 


First, the transport in the active material is considered. At any point of the 
electrode, lithium is intercalated into the active material. The inserted lithium is 
assumed to obey a Fickian type diffusion transport process which can be written 
as 


дс, 
дї 


where c, is the lithium concentration апа D, is the corresponding diffusion 


=V. (Dyc), (4.1) 


coefficient of active material. In addition, the active material secondary particles 
are assumed to be spheres such that spherical coordinates can be employed. 
Moreover, the transport is modeled as spherically symmetric, such that, using 
the transformation operation from Appendix A.1, Equation (4.1) can be rewritten 
in spherically symmetric form yielding 


dc, 109 (xD =) | (42) 


ar y? dy 5 ду 


Note that the dimension у is the radial coordinate inside the spherical particles. 


Conservation of electronic charge in the solid phase 


In the porous electrode, the electronic charge is carried by the solid phase. The 
microscale electronic current density is described using Ohm’s law as 


i, = KEV ps, (4.3) 


where кё" 


son is the electronic conductivity and Ф; is the electronic potential of the 


solid phase. Note that, analogous to Section 2.2, the coordinate & refers to the 
microscale. Due to conservation of charge, the continuity equation reads 


Ve i, =Ve- (кету: ф,) 0. (4.4) 
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From Section 2.2, using Equation (2.24), the macroscopic form of Equation (4.4) 
is 
у. (куф, - as. =0, (4.5) 


where, comparing to Equation (2.24), k Ken is the effective electronic 


ое = 
conductivity, Dg, = Ф, is the macroscopic electronic potential and agg = ase is 
the specific surface area of the solid phase. Electronic charge is produced due to 
electrochemical reactions at the interface between active material and electrolyte. 
This is accounted for by the reaction term bag = j,.F which represents the 
exchange current density induced by lithium flux j,, at the interface between 


solid and electrolyte phase, where F is the faraday constant. 


Rearranging Equation (4.5) leads to 
у. (KRV) easi. (46) 


which is the macroscopic continuity equation representing the electronic charge 
conservation in the solid phase. 


Conservation of ionic charge in the electrolyte phase 


The ionic charge is carried by the electrolyte phase. Using the concentrated 
solution theory on binary electrolytes [47], the ionic current density can be 
written as 


> 


i, = KP" V. 9, — k5" V. Inc,, (4.7) 


ion 


where K7?" is the ionic conductivity, Фф, is the ionic potential and c, is the lithium 


concentration of the electrolyte phase. Moreover, the diffusional conductivity is 


"m VERT | s, + B SOC, ( $ E) | (48) 


Bp F пу тУ+ nco д1пс„ 


where 54+, so and v are stoichiometric coefficients of the cations in the solute, z+ 
is the charge number of the cations, " is the transference number of the cations 
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with respect to the solvent, n is the number of moles of a species and cg is the 
concentration of the solvent. In case of LiPFe, which is the electrolyte salt used 
in this work, the reaction equation is LiPFg == Lit +РЕ [10] and, therefore, 


the parameters are 5; = —1, sọ = 0, n = 1, v; —1,v—2andz, = 1. Note 


that v = v} + v_, where v_ = 1 15 the a stoichiometric coefficient of the anions. 
Equation (4.8) can therefore be simplified as 


up. SKORT aln fa 
ion _ e 0 | EN И 4.9 
Е: F (+ ) ( olnc, ) a 


Due to conservation of charge, the continuity equation reads 
Vid, = Ve (Vege) + Ve (XB Vene) =0. (4.10) 


From Section 2.2, using Equation (2.25) and recognizing that the electrolyte 
represents the B-phase, the macroscopic form of Equation (4.10) is 


М» (=к УФ, ) + у= (KB, V Inz.) + ase j seF =0. (4.11) 


where, comparing to Equation (2.25), the effective conductivity kg eft 18 identified 
as either the effective ionic conductivity коп. or the effective diffusional 


conductivity 


10n 
" AKT T ( i am) (a 1) | (4.12) 
i F alng, 

Moreover, the macroscopic potential pg is either the macroscopic ionic potential 
@, or the macroscopic concentration of lithium c, of the electrolyte phase. 
The specific surface area laß = ase is between the solid and the electrolyte 
phase, where ionic charge is produced due to electrochemical reactions. This 
is accounted for by the reaction term j F, as well. Rearranging Equation (4.11) 
yields the macroscopic continuity equation representing conservation of ionic 
charge in the electrolyte phase as 


у. (котур, – ken V Ine, = аш}. (4.13) 
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Conservation of mass in the electrolyte phase 


Following the above assumed concentrated solution theory of binary 
electrolytes [47], the flux density of cations is 


> 


dinco 19 " 
) Vece + LF cio, (4.14) 


dinc, 


Jı=-vD, (1- 


where c, is the salt concentration of the electrolyte or, in short, electrolyte 
concentration, D, is the diffusion coefficient co is the concentration of the solvent 
and Vo is the velocity of the solvent. The corresponding continuity equation 


reads 
Qc, 3 
rix 
dinco id 
=V.-[v;D, (1 V et V 
5 | Al im) EHE en| (4.15) 
Ve- (i,t? 
dincg é ( e e) 
= У; 0. ( = Te) Ve | = "ES -Ve C4 V0 
Using electroneutrality conditions, i.e. c, = с+/у+ =c_/V_ or с+ — c,V4, 


Equation (4.15) can be rearranged as 


дс dInco Ve: (10) 

v — = уу: |р, {1 У — ——_——~ — ViVe- (c,Vo): 

E uc: 4 im) 56е uF +Уа (с) 
(4.16) 

Dividing by v+ on both sides brings 
dc, dinco iVe au 

=V;-|D.{ 1- V — = Ve- (c, 4.17 
дї & | il im) sed aVF ё (с.о), ( ) 
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0 

where the identity Vz - mu =з) +i,Ve -t9 was used. The first 
part of Equation (4.17) is due to the diffusion process, where D, is the diffusion 
coefficient and the driving force is the concentration gradient Vzc,. In the second 
part, the migration part, is governed by the electric current density i. Finally, the 
third part is due to the velocity vo of the solvent and is called convection part. It 
is usually assumed that dIn co/ dInc, ~ 0 and vo = 0 [36, 46]. Additionally, i? is 
treated as concentration independent, which was observed by experiments [130] 
and which makes the gradient term vanish. However, in [131], it was indicated 
that this may not be correct. Finally, the microscopic conservation of mass of the 
electrolyte phase can be reduced to 


дс 

zt = Ve (Б„Уес,). (4.18) 
From Section 2.2, using Equation (2.25), the macroscopic form of 
Equation (4.18) is 


дс _ м 
9437 =V. СУ) ы) а, (4.19) 


where, comparing to Equation (2.25), фр = ф, is the volume fraction, pg = с, 
is the macroscopic concentration, kg eff =D, eff is the effective diffusivity and 
Aap = Ase is the specific surface area of the electrolyte phase. Recall that the 
reaction term je represents the lithium flux density at the interface between solid 
and electrolyte phase. Note that it is expected that the anions are not taking part 
in the electrochemical reaction which is accounted for by the term (1 — ©). Thus, 
Equation (4.19) represents the macroscopic continuity equation describing mass 
conservation of the electrolyte phase. 
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Reaction kinetics 
The electrochemical reaction at the interface of electrolyte and active material 


is described by a Butler-Volmer type equation [10, 47, 132] which, in case of a 
lithium-ion cell, can be written as 


= P" io (1-a)F_ aF _ 
Jse = Ce) = т є exp RT n exp RT” 


= —koc, © (Coase = С өш) Cs surf 
oe (l-a)F_ = OF _ 
PM 7) ер 2—17) 0, 
(4.20) 


where ig is the exchange current density, ko is the effective reaction rate 


constant [10], @ or (1 — о) are symmetry factors representing a favouring 
of cathodic or anodic reaction, respectively. Typically, @ is set to 0.5 [47]. 
Additionally, c, max is the maximum concentration in the solid phase and c, sur is 
the concentration at the particle surface. Here, 


UE (Ф, = 9.) = Ега (Cs surf) (4.21) 


is the overpotential indicating if the potential difference between solid and 
electrolyte is above or below the equilibrium potential Elan) The latter 
of which is a function of the concentration at the surface of the particle. 


4.1.2 Classical half-cell model 


In this Section, the half-cell model is presented as it was implemented in this 
work using COMSOL Multiphysics [133]. Using the software’s equation-based 
modeling module, partial differential equations (PDEs) can be entered directly. 
In Figure 4.2, the implemented PDEs are summarized, which are the previously 
presented macroscopic continuity equations. 
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The transport is considered through thickness direction of the cell only. Based on 
this assumption, transport in the lateral directions is neglected. Therefore, only 
one-dimensional transport is considered. Figure 4.2 shows two model levels. 
Each of the levels contain one-dimensional domains. These are, on the one 
hand, the separator and positive electrode domain, i.e. cathode, on the cell level. 
The separator and cathode, i.e. positive electrode, domain is denoted by the 
superscript "sep" and "pos", respectively. Note that the conservation of charge 
in the solid phase is omitted in the separator part, since there is no solid phase 
present. On the other hand, the one-dimensional particle domain can be found on 
the particle level. 


On the top level, the cell level, the macroscopic continuity equations, i.e. 
Equations (4.6), (4.13) and (4.19), for the respective cell quantities of interest 
are defined. The respective cell properties of interest are the electronic potential 
in the solid phase @,, the ionic potential in the electrolyte phase @, and the 
electrolyte phase concentration c,. On the particle level, the mass conservation 
in the solid phase, i.e. Equation (4.2), is present. The cell property of interest is 
the concentration of lithium in the solid phase c,. 


Boundary conditions 


In addition to the partial differential equations, the boundary conditions (BCs) are 
also shown in Figure 4.2. In the special case presented here, the anode is modeled 
as a lithium metal foil. This assumption leads to different boundary conditions 
compared to the full-cell model, see [36, 42]. Furthermore, galvanostatic 
discharge conditions are considered, thus, a constant electric current is applied. 
In case of galvanostatic charging processes, the sign of the applied current has 
to be reversed. As was mentioned before, transport in only one direction, i.e. x- 
direction in Figure 4.2, is considered such that volume integrals become simpler, 
see Appendix A.2. 
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Electronic charge in the solid phase 


As a first step, the total electronic current in the solid phase is calculated. 
Therefore, Equation (4.6) is integrated over the cathode volume yielding 


gell ET V к PSVo dx = Ace! pm E Жах 
| K, eff Ф, SE sep Ase] se 


Lsep 

0 

А Ki УФ, CEE — де ae Ф, (L*P)) = 

tot 

Ac Pd азе), dx, 
(4.22) 
where 

кт "уф, (1%) = 0) (4.23) 


accounts for the fact that the electronic current density must be zero at the 
separator-cathode interface. Moreover, the electronic current density at the 
cathode-current-collector interface is identified as 


tot = 
EVD L) = [^ asi, Fir. (4.24) 
SHEE x 


sep 


i, ( ztot) 


Due to galvanostatic boundary conditions, the electronic current density at the 
cathode-current-collector interface is set as 1,(//9©) = ip. Note that the total 
cell length is the sum of the separator and the positive electrode length, i.e. 


ЈА — [sep + 1Р8. 
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Concentration in the solid phase 
Next, the volume average change in concentration of lithium in the solid phase is 


calculated. Therefore, Equation (4.2) is integrated over its spherical domain and 
averaged over its volume, see Appendix A.2, which leads to 


3 "sec 596; 3 | д ( 5 дс ) 
x At” n3 Oy OS ay J 


Гес f ду sec : oy | 


where r,.. 18 the radius of the spherical active material particles. Note that 


9c, ес 5 
TE 


D, д у $е 


(4.26) 
is set to account for the electrochemical reaction boundary condition at the 
surfaces of the particles [40, 42], where j,, is the Butler-Volmer type flux density 
from Equation (4.20). In other words, by inserting the boundary condition 
from Equation (4.26) into Equation (4.25), the volume-averaged change in 
concentration can be expressed as 


1 "ес дс, = 
zf y э, Ф = е (4.27) 
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lonic current and potential in the electrolyte phase 


The total ionic current in the electrolyte phase is calculated by integrating and 
volume-averaging Equation (4.13) over the cell volume as 


15°Р 
де [^ v (Iv, Kar Vinc.) de 
tot 
cel J v.(- ктеу, — koe PV Ing, ) dx = (4.28) 
Lsep 


tot 


EC f 2 Ase} oF Ax. 


The left-hand side is 
ll bs 
к!" se] ion, se] 
А® n Ve (- K, eff РУФ, – Kp eff Рус, 2 dx+ 
tot 


L 
cell yon ‚pos ion,pos E E 
А is NS (- K, etf УФ, — Kp eff Vinz,) dx = 


cell gen Sep se ien;sep = (те 
A (- K, efr У r = Kp eff Vine, (L °))- 


A (—кушг°?Уф,(0) — кру Vinz,(0)) + 
0 


(429) 
ae (eve uy е Ie )- 
A (KV, "Vinc (L99)) = 
— дее! (- кү "V9, (0) 2 Kr Vn c.(0)) 
where 
кіну (L) — PSY Ing, (L) = 0 (4.30) 
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accounts for the fact that no ionic current is allowed at the cathode-current- 
collector interface. Also, continuity of ionic current at the separator-cathode 
interface is represented by 


___lon,sep sep) __ ,.ion,sep = sep) _ 
Ko eff V@,(L*?) Kp eff Vine, (LP) = 


ion, pos sep ion,pos = (TSep (4.31) 
— Ke eff V@.(L hr Kp, eff VInc,(L ). 
Finally, using the results of the left-hand side simplifications, leads to 
[tot 
lon,sep ion,sep _ E 
"CH ve 0) -крш Уш) = f. aee Pes qa 


where the total ionic current density at the anode-separator interface is expressed 


as 
[tot 


10) = | aci, Far. (4.33) 
Lsep 

Due to galvanostatic conditions, and by comparing Equation (4.33) to 
Equation (4.24), the ionic current density at the anode-separator interface is 


Concerning the ionic potential, usually, the focus is on the potential difference, 
thus, Ф, is arbitrarily set to zero at the anode side [40, 42], which is expressed 
via 

ф.(0) 20. (4.34) 
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Concentration in the electrolyte phase 


The total flux in the electrolyte phase is achieved by integrating Equation (4.19) 
and averaging over the cell volume, which results in 


15°Р tot 
a А sep IC dx 4- op le 2) 


Lot 3 im 
І T sep w 
Som f v (рус, ax 
[tot _ 
ur v (orav) оаа) (4.35) 


1 ET pP nr pos w= 
Е Lot [ P (р ett Се )&+ тер У. (Dive. я 


1 ГАС 
n ase (1 — 19) js dx. 


Lt | sep 


The divergence terms give 


15°Р Е 
р OS 
), у. (р ыг Се )&+/ v ‚(ре Уе, ) ах 
0 
= DENE?) — Doe Ve, (0) + Р VEAL) — Рр Vet”) 
= -DE УС, (0) D 
(4.36) 
where 
Di eV Ce (LP) = Dorit’ Ce (LSP) and ЮРУ Ce (129%) =0 (4.37) 


reflects continuity of mass flux between separator and cathode region and implies 
that no mass flux is allowed at the cathode-current-collector side, respectively. 
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Using Equation (4.36) in Equation (4.35) yields 


1 LSEP di dc T ‚ztot 2 
= pe dy pos Се ах 
ЈА" | 0 е ot хер Фф; 


-L (р ve e(0)) — ET , ^ (1—19)7 dx 
[iot (р e,eff ) [iet jn se +)Jse 

(4.38) 
Equation (4.38) shows that the volume-averaged temporal change of the 
electrolyte concentration is due to diffusion at the anode-separator interface and 
the electrochemical reactions of cations at the electrolyte-solid interface. As 
was assumed earlier, the applied electroneutrality condition guarantees the same 
amount of anions and cations in a region. Globally, this is expressed as a constant 
average electrolyte concentration, i.e. a vanishing temporal change. Therefore, 
the left-hand side of Equation (4.38) is set to zero yielding 


tot u 
рн, (0) = / „ы e Hise d. (4.39) 


In other words, the flux of cations due to electrochemical reactions at the surface 
area of electrolyte and active material is compensated by migration at the anode- 
separator interface [38]. Comparing Equation (4.33) to Equation (4.39), the right- 
hand side of the latter can be rewritten in terms of the applied current density as 


tot int (1 _ t?) [et z ls А 
у Ase(l =t} )isedx = | Ase] so-F dx = = (1 і), 
—— 


(4.40) 


i (0) =lapp 


where, following the assumptions made in Section 2.2, the specific surface area 
a,. and, additionally, i are assumed constant and, therefore, can be extracted 
from the integral. Thus, the boundary condition for the flux in the electrolyte 
phase at the anode-separator interface can be written as 


i 
рус, (0) = 1-0). (441) 
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The boundary condition in Equation (4.41) is commonly used to replace the 
anode in a half-cell model [40, 42, 132]. However, disassembling the continuity 
equations for the electrolyte phase, as was done above, shows the implications 
that come along with it. Such as, the volume average lithium concentration in the 
electrolyte can be assumed to stay constant over time. 


Cell quantities 


The half-cell model, as presented above, is used to evaluate cell performance. 
To this end, typically, cell quantities like state of charge, specific capacity, cell 
voltage and alike are computed and discussed. 


Specific capacity 


The total capacity, or ampere-hour capacity, or electric capacity, of the cell can 
refer to either the negative or positive electrode. It is calculated via the initial and 
maximum lithium concentration in the solid phase of the respective electrode. In 
case of the half-cell model used in this work, the total capacity of the positive 
electrode is given as 


Ош = VP% (ar E am 
cell jos (4.42) 

=А Г, Фах u Cmts 
where УР" and @, is the volume fraction of the solid phase of the positive 
electrode, respectively. A“! and LP® is the volume, cross section area and length 


of the positive electrode, respectively. c and c, max is the initial and maximum 


sinit s,max 


allowable concentration of the active material, respectively. The specific capacity 
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is commonly used in the battery community. For this purpose, the capacity is 
related to the mass of the active material, i.e. mam, which brings 


Il 5 
Q = Ош Ж a LP (e s max m^ F 
spec — MAM = Acel L pos jj. OAM (4 43) 
= (Cs miax E Cia) 
OAM | 


where одм is the density of the underlying active material. 


Cell voltage 


The cell voltage is defined as the voltage drop between positive and negative 
electrode-current-collector interface. Since the anode is represented as a lithium 
metal foil, the anode potential is assumed to be zero. Thus, the cell voltage is 
equal to the cathode-current-collector interface as 


сеп = Ф,(1/®*). (4.44) 


State of charge and depth of discharge 


Typically, the cell voltage is drawn versus the so-called state of charge (SOC) or 
depth of discharge (DOD) of the cell. SOC is defined as the releasable capacity 
relative to its maximum allowable, whereas DOD is defined as the released 
capacity relative to its maximum allowable [134, 135]. In case of a full-cell setup, 
the state of charge or depth of discharge can be related to the capacities of either 
the positive or the negative electrode. As for the half-cell model in this work, 
the state parameters are defined using the positive electrode capacity. Since the 
electrode capacity is directly related to the amount of lithium stored inside the 
active material, the DOD can be calculated as 
Cs avg (t) E Cs init 


DOD(t) = Ha siit] (4.45) 


Cs max — Су init 
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where C, avg (t) is the cell's volume-averaged lithium concentration of the solid 
phase at time f. During discharge, as the lithium concentration increases in the 
positive electrode, the value of Equation (4.45) increases from 0 to 1, representing 
an initial concentration state and a maximum concentration state where, ideally, 
all sites of the active material crystal structure are filled with lithium [10]. In 
other words, the cell's depth of discharge increases while, by the same time, its 
state of charge decreases from 1 to 0. Therefore, the cell’s SOC and DOD are 
related via 

SOC —1—DOD. (4.46) 


Calculating the nC current 


In the battery community, it is very common to use the concept of C rate. The 
C rate represents an electrical current, where a C rate of 1 corresponds to a 
current that is needed to charge or discharge a LIB within one hour. Higher 
C rates stand for higher currents whereas lower C rates mean lower currents. In 
reality, the value of the C rate is highly dependent on the structure, geometry, 
electrochemistry, etc., such that, usually, experiments are needed, in order to find 
the electric current representing a C rate of 1. 


In the cell model, however, the 1C rate representing electric current can be 
calculated directly from the given structural parameters. It was shown that the 
state of charge can be calculated by knowing the state of lithium concentration 
of the electrode. In case of the presented half-cell model, the volume-averaged 
temporal change of lithium concentration of the positive electrode is generally 


дс, аур 1 Let 3 "sec 2 дс, 
savg _ 25 dydx. 4.47 
dr Ір In E Tor en 


Recall that, from Equation (4.25), the temporal change of lithium of a single 


calculated as 


particle is known as 
3 "с „IC ui 
ee Ja ЧУ = e (4.48) 


sec 0 sec 
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Volume-averaged integration of Equation (4.48) over the electrode volume leads 
1 tot 3 Tsec dc 1 [tot 3 E 
n sf y? —3 dydx = a= | se (4.49) 
LPOS Jrsp Fec? JO ot LPS Jrsep Tee 
Inserting Equation (4.49) into Equation (4.47) brings 


дс, 1 tot 3 _ 
= | NET (4.50) 
Lsep Гес 


to 


Additionally, from Equation (4.24) and the boundary condition i i,(L™), the 


app — 
applied current density is represented as 


tot 
"INEST (4.51) 
Lsep 


Under the assumption that the interfacial area a,. and the active material particle 
radii г are constant, they can be extracted form the integrals of Equations (4.50) 


and (4.51), converting them into 


Km 1 3 fe”. 
AME x | Dk (4.52) 
sec 
and 
Lot 
| j,, dx. (4.53) 
хер 


Finally, comparing Equations (4.52) and (4.53), brings 


ee c 1 


ш 3 Im (4.54) 
дї DPSS ты. da 


From Equation (4.54), it can be seen that the volume-averaged temporal change 
of lithium concentration of the positive electrode can be computed by the applied 


current density i,,,, and the structural properties LP95, r,.. and ase. Note that those 


app 
properties are taken as constants in this work, which also makes Equation (4.54) 


constant. 
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In order to calculate С) at time т, Equation (4.54) is multiplied by t and ће 


initial concentration c is added such that 


sinit 


дс 
Cs ave (t) = LO ETE Cs init 
; 4.55 
ae (4.55) 


s= et +C.. 
snit 
Lpos Гес ауе” $,1nil 


Inserting Equation (4.55) in Equation (4.45), the depth of discharge at time f is 


1-3. ар. 
Dodge memes. (4.56) 
Cs max — Cs init 


A fully discharged cell is represented by setting DOD (t) = 1 within one hour 
t = 3600s in Equation (4.56). Rearranging yields the 1 C-current as 


ac LP98y Ase F 


lapp = (Сая Cs init 3 73600 (4.57) 
Generally, the nC-current is 

a LPS p ec Use F 

iapp = n(Cs max Cs init 3 33500 (4.58) 


Note that the above derivations do only apply for cases where the specific 
surface area can be calculated using ase = ©, which implies the assumption 
that the active material particles are of spherical shape and detached from each 
other [136]. Morphological variations from this formula necessitate special 
treatment of the above presented derivations, which shall not be considered in 
the present work. 
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Figure 4.2: Mathematical model of the half-cell setup. 
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4.2 Hierarchically structured half-cell 


As mentioned in Section 1, the performance of LIBs can be increased by so- 
called hierarchically structured cathodes where, by postprocessing the initial 
active material, the porosity of secondary particles is increased. As a result, the 
rate capability and cycle stability can be improved [26-29]. In order to simulate 
such structures, in the following, the previously presented half-cell model is 
extended for hierarchically structured cathodes. Consequently, the model is 
called hierarchically structured half-cell model [137]. 


In Figure 4.3a, the sketch of the hierarchically structured half-cell is presented. 
Similar to the classical half-cell setup from Section 4.1, on the left-hand side, the 
anode is a lithium metal foil and, on the right-hand side, the positive electrode 
is a porous composite structure which is composed of the active material and 
the filler material. The latter of which is a carbon-black-binder-mixture. Also, 
a separator is placed between anode and positive electrode. In the hierarchically 
structured electrode case, however, the active material secondary particles— 
which are built up by primary particles—show a distinct porosity. Both the 
separator and composite electrode structure—including the secondary particle 
pores—are filled with liquid electrolyte. 


During charge or discharge, electrochemical reactions take place only at the 
interfacial area between primary particles and electrolyte inside the secondary 
particles. Due to the porous secondary particles, an additional transport region 
is introduced into the system. Therefore, as an extension of the classical half- 
cell model [36, 40-42], the hierarchically structured half-cell model is divided 
into three levels. The cell, the secondary particle and the primary particle level. 
Following the classical half-cell model, the porous electrode theory [35, 47] does 
also apply to the hierarchically structured half-cell model. 
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Anode  Separator Cathode 


222 


Primary particle 
level 


Figure 4.3: Hierarchically structured half-cell setup of lithium-ion batteries. a) Sandwich structure of 
a lithium-ion battery cell. The anode is a solid lithium metal and the cathode has porous 
active materials, i.e. secondary particles. b) The porous structure of both the cathode and 
the secondary particles are homogenized. 


On cell level, transport is carried by the electrolyte and the solid phase, where, 
similar to the half-cell model, the solid phase accounts for the active material and 
the filler material. On secondary particle level, all transport is via the electrolyte 
or the primary particles, where the latter of which is identified as the solid phase 
on this level. Finally, the intercalation and transport process of lithium inside the 
solid phase is modeled on the primary particle level. 


A model for porous secondary particles was also presented in [53]. In 
the present work, however, the secondary particle and electrolyte share the 
same specific surface area for exchange flux. Additionally, to avoid double 
counting of exchange flux density, the electrochemical reactions are only between 
primary particles and the surrounding electrolyte within the secondary particle 
and electrochemical reactions between secondary particle and electrolyte are 
neglected. Moreover, transport of lithium inside the secondary particles is only 
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via the internal electrolyte phase and is intercalated via electrochemical reactions 
into the additionally introduced primary particle level. 


4.2.1 Macroscale equations 


In this section, the mathematical formulation of the hierarchically structured 
half-cell model is derived. The electronic and ionic charge transport, as well 
as cationic flux has to be accounted for on both, the cell and the secondary 
particle level. Therefore, in total, six continuity equations have to be defined and 
solved. Similar to Section 4.1.1, microscopic transport equations are defined for 
the respective phases and volume average theorems are used to convert them to 
macroscopic forms. Additionally, on primary particle level, the lithium transport 
is modeled by a Fickian type diffusion. A summarized version of the presented 
PDEs can be found in Figure 4.4 and 4.5. In the following, the superscript "sec" 
and "prim" refers to the considered level, i.e. secondary particle or primary 
particle level. If no superscript is given, the properties refer to the cell level. 


Primary particle level 


On the primary particle level, lithium moves due to diffusion. Spherical 
symmetry is assumed for the primary particles. 


Conservation of mass in the solid phase 


The transport process inside the primary particles is taken as one-dimensional 
Fickian diffusion by 


дсрі" 19 " дс? 
ðt 292 (: D; dz | | mom 
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Here, the concentration of lithium is ст and the diffusion coefficient in the 
solid phase is D,. The dimension z is introduced as the radial coordinate inside 


the primary particle. 


Secondary particle level 


On the secondary particle level, the derivation of the macroscale forms of 
the conservation equations is quite similar to the derivations in Section 4.1.1. 
Additionally, and as can be observed by experiments [29], hierarchically 
structured secondary particles are of spherical shape. Therefore, spherical 
symmetry is employed and the derived macroscale equations are converted into 
spherically symmetric forms according to Appendix A.1. 


Conservation of charge in the solid phase 


Similar to Equation (4.6), the macroscopic form of the conservation of charge in 
the solid phase of a secondary particle is 


V. (- gu ves = TSF, (4.60) 


se 


eon,sec 


where K; eff is the effective electronic conductivity of the solid phase, Фф is 


the macroscopic electronic potential and а% is the specific interfacial area of the 
primary particle surfaces and the electrolyte inside the secondary particles. Note 
that the macroscopic surface flux of the secondary particles de can be calculated 
using the Butler- Volmer relation according to Equation (4.69). 


Converting to the spherical coordinate system, and employing spherical 
symmetry, yields 


1 д oo sec Е: 
sx vee ja (4.61) 
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Here, the dimension y is introduced as the radial coordinate inside the secondary 
particle. 


Conservation of charge in the electrolyte phase 


Analogous to Equation (4.13), the macroscopic form of the conservation of 
charge in the electrolyte phase of a secondary particle is 


ion,secY—sec ion,sec _зес\ _ secsec 
v.(- xu ver xu Vinas) = -ak FF, (4602) 


where E and EI is the effective ionic and diffusional conductivity of the 


electrolyte phase, respectively, 97 is the macroscopic ionic potential and с% is 


the macroscopic electrolyte concentration inside the secondary particles. 


The spherically symmetric version of Equation (4.62) is 


lo (- 2 ,ion,sec OQ. 2 jion,sec дш n ) = sec 7860 7r : (4.63) 


у? dy e,eff ду =Y Kp ett ду Ase Ise 


Conservation of mass in the electrolyte phase 


Finally, on the secondary particle level, the conservation of mass in the electrolyte 
is derived as follows. Adapting the macroscopic form of the conservation of mass 
in the electrolyte phase from Equation (4.19) as 


sec de^ Sec үу—$ес sec 0 sec 
е Or =V. (DEVE; ) — Use (1 = ге , (4.64) 


where D? is the effective diffusion coefficient of the electrolyte phase inside 


the secondary particle. 
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Converting Equation (4.64) to the spherical coordinate system and employing 
spherical symmetry yields 


dc 


19 дс 

sec e 2 psec e as -sec 

—— =; (yD "a eg | 4. 

e Or yay o ee E) -a e - 5s (4.65) 


Cell level 


In order to derive macroscopic continuity equations for the cell level of the 
hierarchically structured half-cell model, volume averaging methods are applied 
according to Section 2.2. 


Charge conservation in the solid phase 


The microscopic continuity equation regarding the conservation of electronic 
charge on the cell level is analogous to Equation (4.4). Using the volume average 
approach by employing Equation (2.24) brings 


у. (куф) = фс, (4.66) 


eon ; : . QA. ET омо 
where ky orp = Куе is the effective electronic conductivity, Pa = Ф, is the 


macroscopic electronic potential. In contrast to the common half-cell model, 
the electrochemical reactions at the interface between the secondary particle 
surfaces and the surrounding electrolyte, i.e. the surface term in Equation (2.24), 
is neglected in this work. Rather, the production of electronic charges on the 
cell level is due to electrochemical reactions from within the secondary particles 
on the surface area of the primary particles. This is expressed as the volume- 
averaged source term being identified as Qaba = Psecls,sec» where фә is the 
volume fraction of secondary particles without internal porosity and ТЕ is 
the volume-averaged production term of electronic charges inside the secondary 
particles. 
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Charge conservation in the electrolyte phase 


The microscopic continuity equation regarding the conservation of ionic charge 
is identical to Equation (4.10). Using the volume-average approach, by applying 
Equation (2.25), yields 


У-(—к Уф, — kits V INZ,) = фес. (4.67) 


Comparing to Equation (2.25), the effective conductivity, kg tes is either the 
effective ionic conductivity, коз, or the effective diffusional conductivity Kp eff 
on the cell level. Moreover, the macroscopic potential, pg, is either the 
macroscopic ionic potential, ®,, ог, c,, the macroscopic concentration of lithium 
of the electrolyte phase. Additionally, the generation of ionic charges is due 
to electrochemical reactions on the surfaces of the primary particles inside the 
secondary particles. Therefore, the volumetric production term on the right-hand 
side of Equation (4.67) is defined by the volume-averaged production term of 
ionic charges from within the secondary particles md Note that the reaction 
term in Equation (2.25) on the surfaces of the secondary particles is neglected. 


Mass conservation in the electrolyte phase 


The derivation of the macroscopic continuity equation regarding the mass 
conservation of the electrolyte phase starts with the microscopic form provided 
by Equation (4.18). Using volume-averageing, by employing Equation (2.25), 
results in 


дс = 
€ V. a ; 4.68 
e дї M (De. Vz.) + seco sec - oe) 


where, comparing to Equation (2.25), фр = 6, is the volume fraction, pg = c, 


is the macroscopic concentration and kg eff = D, is the effective diffusion 
coefficient on cell level. The electrochemical reaction on the surfaces between 
the secondary particles and the electrolyte is neglected. Instead, the production of 
cations Jis is governed by electrochemical reaction from within the secondary 
particles. 
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Reaction kinetics 


Similar to Section 4.1.1, the electrochemical reaction at the interface between 
primary particles and electrolyte is described by a Butler-Volmer type equation 
on secondary particle level. In case of the hierarchically structured half-cell 
model, it reads 


sec __ :(туѕес „ѕес „вес 
J se --j(n Ce ee) 


= Say, ex (b OF sec ex AF sec 
ee MERT 


= —ko (cz) dio (Cs ings m Co adit) Er (Е) б 


‚4 ex (1—0) ех OF sec 
р RT n p RT” , 


where, compared to Equation (4.69), iy is the exchange current density on 


sec 
5,SU 


particle surfaces. Moreover, 77° refers to the overpotential at the surface of a 


secondary particle level and се; is the lithium concentration at the primary 


primary particle inside a secondary particle by 


т =, (Фу = Фе) m Ега (Cua) E (4.70) 


sec 


where the equilibrium potential, Esq (CH surf); is a function of су. 


4.2.2 Hierarchically structured half-cell model 


As can be seen in Figures 4.4 and 4.5, the hierarchically structured half-cell is 
divided into three levels, each of which containing one-dimensional domains. 
On the top level, the cell level, the separator and the positive electrode domain, 
ie. cathode, is modeled. Here, the superscripts "sep" and "pos" indicate the 
respective region. 
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Boundary conditions 


Similar to the half-cell model from Section 4.1.2, the anode is modeled as 
a lithium metal foil. Again, galvanostatic charge or discharge is considered, 
thus the applied current is assumed to be constant. In the following, boundary 
conditions are derived by integrating over domains in the cartesian and spherical 
coordinate systems. The mathematical basics of the integration operations 
performed can be found in Appendix A.2. 


Electronic charge and potential in the solid phase 


As a first step, concerning the right-hand side of Equation (4.66), the volume- 
averaged source term Īss is evaluated by volume-averaged integration of 
Equation (4.61), yielding 


3 бес д 2 „.eon,sec 9%% ) 3 sec 2a <sec 
RE AA K d == vj S d 
eh С ИЕ ж 


ѕес 


2 ,.eon,sec ops ( Rec) _ 


— K EE IPEA u 
ru) "sec Ks eff ay (4.71) 


=N PETF dy, 


where r,.. 15 the secondary particles radius. Consequently, by recognizing that 
the electronic current density at the surface of the secondary particle is 


IPE uu) 1 fe 

eon,sec sec) __ 24 secsec 

K, eff тү TN y ase Де FW; 
Гес J0 


ду (4.72) 
——— 
surf 
liaec can be expressed as 
т 3 sec 3 Tse 
15 ѕес = "EE = r3 fie Ya. F dy, (4.73) 


sec sec 
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which, in turn, can be used for the right-hand side of Equation (4.66) yielding 


ЕЕ or тес 
Psecls,sec = ye ls, surf — Asels surf > 


"sec (4.74) 


dse 


where, in case of equal-sized and detached spherical secondary particles, ase = 
Е [136] is identified as specific surface area between secondary particle and 
electrolyte. 


As the next step, the electronic current over the cell domain is calculated by 
volume-averaged integration of Equation (4.66) as 


tot tot 


cell Keon,pos cell D 
A | у. (- Ko, ) dx = А1 | 6. nsec dx 

LseP LseP 
0 
cell Koon Pos tot cell eon,pos Se P 
A K eff V9,(L FA (ке s e) — 

tet 3 x 
cell "ee 2a вест —sec 
A бес 3 3 б у Ase E dy dx, 


(4.75) 
where А! is the cross section area of the cell. Note that the result from 
Equation (4.73) was used. Additionally, 


eon,pos sep\ __ 

—K; etf V@,(L*?) =0 (4.76) 
accounts for the fact that electrons are not allowed to enter the separator domain. 
Finally, the electronic current density at the cathode-current-collector side is 


ot 3 Fsec 
: 2a 
Kor V, (17%) Е Tsep Ösec 3 Ü ya pa JF dydx, (4 77) 
m д ѕес t 
i, (Lt) 


where 1, (1/9) is set to be the applied electric current density app: 
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Continuity of the electronic potential between secondary particle and the cell 
level is realized by a Dirichlet boundary condition 


ФУ (тес) = Ф, : (4.78) 


Concentration in the solid phase 


Next, the average change in concentration of lithium in the solid phase is 
considered. Therefore, volume-averaged integration over the spherical primary 
particle level of Equation (4.59) brings 


"prim д pii Tprim д д prim 
: 3 f сг de = i 3 i | ёр, = dz 
0 дї о Oz 92 


r prim r prim 


prim rim 
б. 7 дс; (Тр) 3 02 dc? 
= 3° prim “s 5 Г 
r prim dz та dz 


which describes the volume-averaged temporal change of lithium concentration 


0 (4.79) 


inside a primary particle. Note that 


дс?" Geren) sec (4.80) 


s Qz — se 


accounts for the electrochemical reaction boundary condition at the surfaces of 
the primary particles. Here, pd is the Butler-Volmer type reaction term from 
Equation (4.69). 
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lonic current and potential in the electrolyte phase 


The volume-averaged source term on the right-hand side of Equation (4.67) is 
tackled next. Therefore, Equation (4.63) is volume-averaged over its spherical 
domain, which brings 


3 sec д 2  ion,sec д фе 2 jon,sec oln c. 
Eu ду (- CY Re ett ду YK) ett ду y 


3 "sec 
2 а5 sec 


3 Jo y Ase Ise F dy 


Гес 


3 2, ion,sec ope ( Tus) 2 _.lon,sec д Бс Qus) 
Гес En a = 


js 


е, eff ду "sec D,eff 


ion,sec dln се (0) 


3 (- 02x jon,sec x 
3 К, eff D,eff 
"ес I5 dy 


3 Tsec wor 
-= | aie Fay. 
sec 0 
(4.81) 
Identifying that 
ion,sec 9%% (sl ion,sec д In ce (к) = 
e,eff ду D,eff ду 
———————— 
ip surf (4.82) 
1 "sec 2 secsec 


Р Fd 
ЖЕ 0 Y Ase Ise У 


is the ionic current density at ће surface of ће secondary particle, the volume- 
averaged ionic current density is 


т 3 sec 
buo wor] yas i Fay, (4.83) 


sec Vee? 
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where А is the ionic current density at the surfaces of the secondary particles. 
The right-hand side of Equation (4.67) is rewritten as 


= ae sec 
Psecls,sec = bec le, surf — 4sele surf 


"sec И (4.84) 
beas 
dse 
The total ionic current, i.e. volume-averaged integration of Equation (4.67), is 
LSEP 


Av [^ v (Cx V,- oe Vine, dx 


tot 
act (VE, Kon" VInz,) dx = (485) 
L 


sep 
Let 
cell D» 
A Psecle,sec dx. 


Lsep 
The left-hand side is 
SEP 
Av [^ v (кетр, Rn Ving) de 


tot 
cell yon ‚pos ion,pos » = 
a he ve (- Ko eff Vo. Kp, eff Vinz, ) dx = 


cell K’ Sepa 
A (- K, ett Y 


Ace! (- KEY 0) — KPY Ine 0 ) 4 
K, eff Ф,( ) D,eff 0 a: ) (4.86) 


cell кї°їР%% 77 tot 
—A (x e,eff V@(L 


cell „Klon,pos = 
A (- K eff У 2 


А (— i e vig, (0) — rper VInz,(0)) , 
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where 
KEEP VOL") = кр” аг, (09) = 0 aen 


accounts for the fact that no ionic current is allowed to enter or leave the cathode- 
current-collector interface. Also, continuity of ionic current at the separator- 
cathode interface is represented by 


__lon,sepyy> (TSep\ _ „ON,sep = (Tsep\ _ 
Ke efr V@.(L ) Kp eff Ving, (І. ) = 


ion,pos yy- sep ion,pos усы sep (4.88) 
— Ke eff V@.(L ) — Kp eff Ус, (L ) . 
Using Equation (4.83), the right-hand side of Equation (4.85) becomes 
tot tot 3 r, 
I Se —sec 
Awl i Ösecle,sec dx = А! ГА Ф.-с 3 уат Fdydx. (4.89) 
Lsep хер Tec 0 
Finally, the total ionic current density at the anode-separator interface is 
- (CQ v9.(0) - kgs" VInz, (0) = 
—————— 
ї,(0) (4.90) 
г" 3 "sec P. -sec 
- f. без] y ase DM F dydx, 
LSP Гес 0 


where, comparing to Equation (4.77), 1,(0) is the applied current density top: 


Note, the continuity of the ionic potential between secondary particle and the 
cell level is realized by a Dirichlet boundary condition 


Фе (гс) = Ф, : (4.91) 
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Concentration in the electrolyte phase 


As a final step, the volume-averaged source term on the right-hand side of 
Equation (4.68) is dealt with by volume-averaged integration of Equation (4.65), 
which results in 


3 [ry sec gc dy = 3 I д у2 0% dee 
ra Jo e gt fa Jo ду е ду 


вес 
жуа (1—55)7,, dy 


NETUS 
"E sec ^^ e,eff oy 
= 0 
= 3 2 psec де 
m eff ду 
3 ес 2 sec 0 sec 
m 3 Y Ase (1 ei dy. 
Гес” 70 
(4.92) 
The flux density at the secondary particle surface, a is identified as 
sec Be ge) — 1 т 2 „sec 907% а 
eet ду Treo’ Jo 7° t 
> 
Te sur (4.93) 
1 see 2 вес 0 sec 
Е x У бе d se d 
thus 
B 3 sec 3 "sec 5 Ic“ 
Je,sec Z — Je, surf = = 3 f e 5, ЧУ 
sec sec (4.94) 
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where j,.,,¢ is the mass flux density at the secondary particle surfaces, which 
renders the corresponding term in Equation (4.68) as 


= sec 
PsecJesec = du Tac js E = AseJe,surf - (4.95) 
—/ 


dse 


Applying volume-averaged integration on Equation (4.68) and using the result 
from Equation (4.94) leads to 


хер = tot 
i sep дс, didi pos Се дс e dx 
List 0 дї Lsep дї 
xm Lot 
р OS ў 
= Е | Y (р е, eV Се ) dx T e Ys (р РУ Се ) T Psec Je sec IX 
Е a р“ v tot m рр ux ae 
Е a D, eff Y Ce HW Dep р. eff Y Се 
1 ie деѕес 
+ Lot (t. фес c yoe ar dy 
3 


LP Гес 


Гес 
= a ya seni EU) a) 
sec 
1 LSP бер Ж tot 25 _ 
=a Í у. (DR VE) dr + / V (DEVE) d 

1 i 3 ес 2 „sec ger 

EET » —— dydx 

sec "c / уф 2, D 


Ltt J rsep 


1 [tot 3 r 
6 24 sec sec 
—— — 1-1 d 
[ot Í, фес Fac? Ll» I РУЛ, se GY 
(4.96) 
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The divergence terms read 


h 


15°Р tot 


у. (De, Ve.) а n (Dis ve, dx 


0 
= DENE?) — ОзУ, (0) + Doce Vest) — Dg es) 
= РУС, (0), 
(4.97) 
where 


Dm VC (LP) = Dim Ve. (L^) and DV Ce (L'™*)=0 (4.98) 


reflects continuity of mass flux between separator and cathode region and implies 
that no mass flux is allowed at the cathode-current-collector side. 


Using Equation (4.97) in Equation (4.96) yields 


1 LSP dc tot dc 1 
Lit | 0 er дї dan Lsep 3r ; d Let (р eat V Ce (0)) 


1 Ке 3 "sec 2 „sec I 
LACE. жон sec Се dydy 4.99 
List p sec Гес? / У Pe ot y ( ) 


1 re ó 3 Ts xen (1 gy -— duum 
281 Erz а5%(1] — 
Lt \ Jese "98р Jo Y ау y 
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Finally, the volume-averaged temporal change of concentration in the electrolyte 
phase is 


1 LSEP де, Lot әс, 


Lsep 
b f 3 [бе „ „дб 

+ =, pen ѕес "e d dx = 
n 7 / = Prec 3 " уф y V 


as 1 (р р? ус .(0)) = d UM =; [Pass ATE war). 
Lot еї Ltot | /тзер '5©© ЖЕ 0 se se 

(4.100) 
Equation (4.100) shows that the volume-averaged temporal change of electrolyte 
concentration is due to diffusion at the anode-separator interface on cell level and 
electrochemical reactions of cations at the electrolyte-primary-particle interface 
on secondary particle level. As was assumed earlier, the applied electroneutrality 
condition enforces the same amount of anions and cations. On cell and secondary 
particle level, this is expressed as a constant volume-averaged concentration in 
the electrolyte, i.e. a vanishing temporal change. Therefore, the left-hand side of 
Equation (4.100) is set to zero yielding 


[tot 
рес, (0) = Prec 3 EU yas vie (1—-00)j 7 dydx. (4.101) 


Lsep 


By comparing Equations (4.90) and (4.101), the latter of which is rewritten as 
[tot 3 Tse 
sec 
/ 6.5 |“ Ya -0) dydx — 
L Гес 0 


ѕер 
(1—20) = 2 sec lag (1 — 9) 
a n sec "er i ya a j se Fdydx 73 ’ 


C AGED VN эъ. ш, 
i,(0)— lapp 


(4.102) 
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where a is assumed to be constant. Using the boundary condition from above, 
i.e. setting the ionic current density at the anode-separator interface equal to the 
applied current density, yields 


i 
-Deen Ve.(0) = 2201—08). 


Interestingly, Equation (4.103) is the same as for the half-cell model 
[40, 42, 132]. In order to satisfy continuity of concentration, the concentration 
between secondary particle surface and the cell level is prescribed as Dirichlet 
boundary condition 

er =e. (4.104) 


Cell quantities 


Analogous to the half-cell model from Section 4.1.2, the evaluation of the cell 
performance is done by computing the appropriate cell quantities. Note that the 
cell voltage and state of charge or depth of discharge is computed similar to the 
half-cell model from Section 4.1.2 and are therefore omitted here. 


Specific capacity 
The specific capacity is calculated by first evaluating the total capacity as 


Ош = X [Cosi EE Cs init) F 


cell y pos sec (4.105) 
=A" L 0,0; (Санан eu sacs 


where VP?* is the volume of the solid phase inside the positive electrode, and 


€, init and С; max 1S the initial and maximum allowable concentration of the active 
material, respectively. The former of which, can be calculated using the positive 
electrode's cross section area A“®!, length LP’, volume fraction of secondary 


particles ф, and volume fraction of the secondary particles’ solid phase @°°°. 
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Finally, the specific capacity is 


1l S 5 
Q A Ооо = A [Р°% ф, фу (are Е Cyin) F 
spec = = 


MAM Acel g pos Фф, psec OAM 
F 


(4.106) 


Ç 


(c5 fax > 


= j 


QAM 


where oam is the density of the active material. 
Calculating the nC current 
In case of the presented hierarchically structured half-cell model, the volume- 


averaged temporal change of lithium concentration in the positive electrode is 
calculated as 


де, 8 Xe [эвде 
tele ah "ah ar eh. (4107 
дї LPS }үзер КЗ JO Paus 20 дї 
From Equations (4.79) and (4.80), it is known that 
3 Tprim дате 3 — 
De =, (4.108) 
Porim 0 dt Porim 


Importing Equation (4.108) into Equation (4.107) yields 


дс, 1 f£ 3 fe, 3 д 
авт i [ FS дуй. (4.109) 


| 3 
SP sec prim 


Moreover, from Equation (4.77), it was found that 


п“ 3 "sec 2 —sec i 
sec 3 n Уаз Ise F dydx = і, (4.110) 


LseP 
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if galvanostatic boundary conditions apply and i,(L™) is set to a constant applied 


sec 


electric current density i Under the assumption that sec» prim and ase are 


арр" 
constant properties, Equations (4.109) and (4.110) аге rewritten as 


1 Li 3 "sec 2 3 sec 1 3 г" 3 "sec 2-sec 
el, 3 y r J se фух = pas A sf У Ise dy dx 


SP Гес 0 prim prim LP Tee 
9 
ITs avg/ дї 
(4.111) 
апа 
г" 2 sec = 3 "ес 27sec 
Р бес 3 x m ya ae j, se Fdydx = PsecAse F 4 3 y J se dydx. 
LseP LSP Tec 0 
Е (4.112) 
Comparing Equations (4.111) and (4.112) yields 
ITs avg EC 1 3 1 1 1 
дї N g Porim = Psec JA (4.113) 
— арр 
F -Tece ЕЕ f 


where by @,.. = se ec it was assumed that the secondary particles are equal-sized 
and detached [136]. с Cs ave (t) at time t is calculated analogous to Equation (4.55), 
which is inserted into Equation (4.45) such that the depth of discharge is 
calculated as c, avg (t) at time t is calculated by multiplying Equation (4.113) by f 
and adding the initial concentration C; init yielding 


дс 
Cs avg (t) — e n ECS init 
i (4.114) 
— pp 3 3 
с NEUE sec EIE Cy init- 


sec se Porim se 
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Inserting Equation (4.114) in Equation (4.45), the depth of discharge at a time t 
is | 
ар 33. 
рор (t) = F Гесйзе "prim se . (4.1 15) 


С с 


s,max _ “sinit 
The 1C-current can be computed by rearranging Equation (4.115) as well as 
setting f = 3600s and DOD (t) = 1 representing a fully discharged cell within 


one hour as 


pos sec 
L Гес Ase l'prim Ase F 


iapp = (Cs max Cs init 9.3600 (4.116) 
In general, the nC-current is 
ee л (4.117) 
app s,max sini 9.3600 


Note that the above derivations do only apply for cases where the specific 


surface area can be calculated using ase = ots and ае = 3e which 
implies the assumption that the active material particles, either secondary or 
primary particles, are of spherical shape and detached from each other [136]. 
Morphological variations from these formulas necessitate special treatment of 
the above presented derivations, which shall not be considered in the present 


work. 


4.3 Validation of the hierarchically structured 
half-cell model 

The following section aims at validating the previously presented hierarchically 

structured half-cell model [137]. To this end, geometry, structure and material 


properties from [30] were imported into the model and the simulation results 
were compared to measurements reported in [30]. 
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4.3.1 Geometry, structure and transport properties 


In [30], both classical and hierarchically structured electrodes were prepared and 
electrochemically characterized. Moreover, the morphology of hierarchically 
structured electrodes was varied in terms of primary particle sizes, inner porosity 
and secondary particle size. Additionally, statistical image analysis based on 
synchrotron tomography was used to investigate the structural properties of the 
electrodes. 


Concerning classical electrodes, dense pristine powder (p-NMC) was used, 
where Li Ni; /5Mn1/3Co; 302 (NMC) was the active material. The hierarchically 
structured electrodes were prepared using nanostructured powder (n-NMC), 
which was obtained by grinding, spray drying and calcinating the p-NMC 
powder. See [30] for a more detailed description of the process. In both cases, 
slurries were produced by adding conductive filler material to the powders. The 
conductive filler material comprised of polyvinylidene difluoride (PVDF) binder, 
carbon black and graphite. Finally, the electrodes were produces by casting 
the slurries onto an aluminium foil. In Tables 4.1 and 4.2, the geometry and 
structure properties of the resulting electrodes are summarized. Note that the 
specific surface area is computed as a = 30 [136], where—depending on model 
level—the respective value is calculated by 


3 pos 3 psec 
As б or au = Д Р (4.118) 


ѕес А prim 


respectively. 
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Figure 4.4: Mathematical model of the hierarchically structured half-cell setup. 
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Figure 4.5: Mathematical model of the hierarchically structured half-cell setup. (continued) 
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Table 4.1: Geometry properties of the classical and hierarchically structured half-cells taken 


from [30]. 
Parameter p-NMC n-NMC-F850 n-NMC-F900 
L* [m] 260-1076 260-1076 260-1076 
ІР [m] 50-1076 76-1076 71-1076 


In the following, the classical electrode setup named p-NMC is compared to 
hierarchically structured electrodes denoted by n-NMC, which are based on the 
same material. Moreover, the secondary particle’s radii are comparable among 
all electrodes. Additionally, the nanostructured secondary particles were treated 
with different calcination temperatures of 850°C and 900°C, which led to larger 
primary particle sizes and lower porosities in case of the higher temperature. In 
order to distinguish between those two, the n-NMC electrodes are denoted by 
n-NMC-F850 and n-NMC-F900. 


The effective transport properties are calculated using the bulk material property 
multiplied by an effective transport parameter Коң which incorporates the 
morphology of the transport paths of the corresponding conducting phase. On 
the one hand, it can be calculated by employing the resistor network method 
for both the solid and pore networks on virtual particle assemblies, as presented 
in Sections 3.1 and 3.1. On the other hand, the statistically derived so- 
called M-factor can be used equivalently which accounts for volume fraction, 
tortuosity and constrictivity of the transport phase [92, 93]. One prominent 
way of calculating effective transport parameters, however, is the Bruggeman 
correlation [71, 138], which is a function of the volume fraction of the transport 
phase and the Bruggeman transport coefficient brugg in the form of ks = Heuss, 
As for the cell level in this investigation, this type of relation shall be used to 
compute the effective transport parameters of the cathode as follows. 
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Table 4.2: Structure properties of the classical and hierarchically structured half-cells taken from [30]. 


Parameter p-NMC n-NMC-F850 n-NMC-F900 


Cell level 
Р Ја) 0.50 0.50 0.50 
е [-] 0.54 0.57 0.58 
P |—] 0.28 0.32 0.30 
Фф [-] 0.18 0.11 0.12 
asec [m] P 2.00 - 10° 2.31- 10° 2.05 10? 


Secondary particle level 


sec —] 0.00 0.46 0.38 
ф$°° [—] 1.00 0.54 0.62 
rc [m] 4.20 : 10-6 4.15: 10° 4.40: 10° 
a mcis à 9.00- 106 7.75 100 


Primary particle level 


Frim [m] - 0.18- 1076 0.24.1076 
a) Assumed. b) Calculated using Equation (4.118). 


In [92], the RN was extensively used to generate a large database of sphere 
assembly scenarios with randomly distributed and slightly overlapping particles 
following a polydisperse size-distribution. By means of statistical methods, 
prediction formulas for the transport parameters of the considered packings were 
developed in terms of the mean contact angle and the standard deviation of 
the particle radii. Notably, it was found that a Bruggeman type correlation 
can be employed to compute effective transport parameters of the pore phase. 
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In contrast to the transport through the solid phase, it seemed that neither the 
considerably large polydispersity nor the degree of the investigated densification 
had an influence on the transport through the pore phase. There, a Bruggeman 
coefficient of 1.342 was found to achieve the best approximation [92]. 
Considering that the underlying degree of polydispersity and densification 
corresponds to the cathode structures in [30], the transport coefficient is therefore 
applied in this investigation. 


In addition to that, as a simple model assumption, the electronic and ionic 
transport carried by the filler phase and electrolyte phase, respectively, is modeled 
by splitting the above-mentioned Bruggeman-based effective transport parameter 
of the total pore phase according to the volume fractions. This way, both phases 
can be regarded as smeared homogeneously over the pore phase. Recall from 
before that the filler phase comprises of binder, carbon black and graphite. In 
this work, the filler phase is equipped with an effective electronic conductivity 
resulting from all three components. 


In this investigation, using the volume fractions фу" of the filler and ф?°% of 
the electrolyte phase from Table 4.2, the volume fraction of the pore phase is 
obtained as 

роз = фроз + PPOs, (4.119) 


which is used to calculate the effective transport parameter of the pore phase as 


{роя qpos 1.342 (4.120) 


eff,pore pore 


In the framework of the cell models presented in this work, electronic and ionic 
transport is via the filler and electrolyte phase, respectively. Using the additive 
split of the pore phase from Equation (4.119), the effective electronic and ionic 
transport properties are assigned according to the volumetric share of the filler 
and electrolyte phase as 


pos pos 
( eon.pos — apos 1.342. | Фф | айа jion.pos — Apos 1.342. | be | 
eff ~~ Fpore pos eff ~~ Fpore pos ’ 


pore 
(4.121) 
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respectively. In case of the separator on the cell in this investigation, its porosity 
was assumed as 0.5 , which is in the range of the typical characteristics of 
commercial separators, and its Bruggerman transport coefficient was taken as 
3.0, which is within the range reported in literature [139]. 


Finally, effective transport parameters of the primary particle and pore networks 
on secondary particle level for the hierarchically structured half-cell model were 
used in form of the M-factors reported in [30], which are based on the real 
morphology. The effective transport parameters are summarized in Table 4.3. 


Table 4.3: Transport properties of classical and hierarchically structured electrodes taken from [30]. 


Parameter p-NMC n-NMC-F850 n-NMC-F900 


Cell level 
ропров ..] a) 0.48 0.50 0.51 
kea Po [-] 9 0.16 0.10 0.11 


Secondary particle level 


kionsec [_] - 0.33 0.10 
ke" [=] - 0.20 0.46 


a) Calculated using Equation (4.121). 


4.3.2 Electrochemistry and material parameters 


From [30], the equilibrium potential curves, i.e. open circuit voltage curves 
(OCV), were taken as the ones measured at the lowest C rate of C/20, see 
Figure 4.6. 
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О  n-NMC-F900 


Cell Voltage [V] 


0 20 40 60 80 100 120 140 160 
Specific Capacity [mAhg-!] 


Figure 4.6: OCV curves taken from [30]. 


The half-cell models, as presented here, need information on the minimum and 
maximum allowable lithium concentration inside the active cathode material. 
Those quantities are computed here as follows. During the discharging 
process, lithium is intercalated into the active material starting from an initial 
depth of discharge Dop™, Arriving at the maximum capacity, the final 
maximum depth of discharge РОР" is reached. The initial and maximum 
depth of discharge is computed by assuming that—due to the stoichiometry 
of LiNi;/3Mn;/3Co;/302— the theoretical maximum capacity of NMC being 
around Ош» = 278mAhg™!, see Appendix B.1. In [30], it was found that 
the reversibly accessible capacity in case of p-NCM is 158mAhg™! and in case 
of n-NCM-F900 is Оз = 161 mAh g !. In this investigation, the same capacity 
of 161 mAhg"! was used for n-NCM-F850. As a model assumption, here the 
cell is discharged until the maximum theoretical capacity of NMC is reached, 
which means РОР" = 1. This is justified by discussing Figure 4.6. It can be 
seen that at a cut-off voltage of 3V the gradients of the discharge curves tend 
towards infinity. This implies that the material’s maximum capacity is almost 
reached. Additionally, due to the lithium metal as anode, it can be assumed that 
the supply of lithium ions is enough for the cathode material to fully lithiate. 
It should be noted that, during the initial charging and discharging processes, 
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the capacity of the NMC material might not be fully exploited because of the 
formation of a lithium-rich phase at the surfaces of the active material particles 
which hinders lithium diffusion [140]. However, there it was also shown that by 
an extended voltage relaxation step at deep discharge, this phase disappears and 
the capacity efficiency can be fully recovered. In turn, theoretically, this means 
that the provided voltage window is enough to fully lithiate the given material, 
i.e. reaching the maximum theoretical capacity. 


Subsequently, the onset of the simulation, i.e. initial depth of discharge DOD", 
is calculated relative to ОРОЮ", By using the reversibly accessible capacity 
from before, the initial depth of discharge is computed as 


oma _ mot 


spec spec 
Qth.max ? (4. 1 22) 
spec 


D opt = 


which is 0.42 and 0.43 in cases of n-NCM-F850/n-NMC-F900 and p-NCM, 
respectively. Furthermore, the maximum and minimum concentration in the solid 
phase is computed by 


th.ma: 


X 
= max spec ONMC 
Cs max = DOD a 


nd 


F (4.123) 


th.max ONMC 
ar init ~Spec 
Cs init = DOD T ’ 


respectively, where the NMC density of onmc= 4770kg m [141] was used, 
which can be calculated according to Appendix B.2. 


In literature, values of NMC lithium diffusion coefficients can vary from 
1:-107!5m?s7! nee [142-146] and the electronic conductivity 
can vary from 2.2: 10-^Sm^! to 5.2-10-°Sm~! [145, 147]. Additionally, 
in [142, 148], electronic conductivity was measured as a function of temperature 
and lithiation state of NMC. Here, the lithiation state is depicted by у = с/с, max 
in LiNi; /3Mn; /3Coj /302, where—theoretically—y = 1 refers to a fully lithiated 
and y = 0 to a fully delithiated state. 
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In Figure 4.7, the electronic conductivity of NMC is plotted over the lithiation 
state at 30°C, where the experimental values from [142] were approximated by 


the fit-function 
1g (Keen o) = —2.988 у!09 — 2.629 y+ 0.624. (4.124) 


It can be observed that the value of conductivity ranges from = 1 to = 
1-107? S m^! for lithiation states of у = 0.25 to у = 1.00, respectively. This 
result by [142] was acknowledged by [143]. There, additionally, the diffusion 
coefficient and reaction rate for NMC was measured. Therefore, in the following, 
material parameters were taken from [143] and [142] because a complete and, 
thus, more likely consistent set of diffusion coefficient, electronic conductivity 
and reaction rate constant for NMC can be expected. 


lg(« Sic) [Sm] 


0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
yin LiyMn;3Ni1/3Co1/502 


Figure 4.7: Electronic conductivity of NMC over lithiation. Experimental values taken from [142] 
and fitted by Equation (4.124). 


The kinetic reaction rate constant is often treated as a fitting parameter [149]. 
However, in case of different electrode structures, the fitting process could yield 
different reaction rate constants for each variation. In contrast to that, the 
approach in this work is to approximate the effective reaction rate constant Ko 
purely based on the underlying material. By taking the exchange current density 
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from experiments on NMC in [143], the exact same material based constant is 
used for every simulation in this work. This way, the comparability between 
different cathode and cell structures can be improved. 


To compute the effective reaction rate constant, first, the Butler-Volmer type 
reaction equation from Equation (4.20) is employed. Based on this equation, 
the reaction rate constant can be calculated as 


io -1-a 1-а a 
TE = —Кос, (С; тах = Cs surf) Cs surf 
д io 1 (4.125) 
0 = 
1-а 1-94 о 
F Ce (Cima = Cs sut) Cs surf 


Secondly, the experimental value of 7.2 Am ? is used for the exchange current 
density io, which is taken from [143]. From Table 4.4, the mean value of 
maximum and minimum lithium concentration in the solid phase for p-NMC of 
35402 mol m? and the mean electrolyte concentration of 1000 mol m^? is used 
for c, surf and c,, respectively. 0 is taken as 0.5. As a result, the effective reaction 
rate constant is recalculated as 


7.2 1 


= z 1.0- 10719 4.126 
F 100005(49477 — 35402)°° 4947705 er 


0 


Therefore, in all simulations in this work, ko is taken as 1.0-107!°mol7!m?> /s. 


The LP30 electrolyte (1 M LiPFg v/v: 1:1 DMC:EC) was considered, where the 
concentration dependent material parameters were taken from [131]. The initial 
concentration is 1000molm7?. In accordance with Section 4.1.1, rae is set toa 
constant value of 0.23, which is the average value of n for the concentration с, 
between the assumed extreme values 0 and 3000 mol m ?. 
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Table 4.4: Material properties. 


n-NMC- n-NMC- 


Parameter p-NMC F850 F900 Reference 
T [K] 298 - 
onmc [Кет] 4770 [141] 
рор" [| 1.00 - 
= Equation 
DODi™ [_] 0.43 0.42 0.42 ic 13) 
E Equation 
C, init [mol m 7] 21357 20823 (4.123) 
3 Equation 
Cs max [mol m 3] 49477 49477 (4.123) 
C, init [mol m?] 1000 [131] 
D, [2571] 1.00. 10-14 [143] 
Keon [Sm-!] 100 [13] 
Ко от 1] Equation (4.124) [142, 143] 
ko [mol-! m?5s=!] 1.0. 10-10 y eee 
a [-] 0.5 [47] 
10 [-] 0.23 [131] 


4.3.3 Results and Discussion 


The above described material, geometry and structural parameters were imported 
into the classical and the hierarchically structured half-cell model from 


134 


4.3 Validation of the hierarchically structured half-cell model 


Sections 4.1 and 4.2, respectively. In the following, discharging currents 
corresponding to different C rates were applied in each model, where the C 
rates ranged from rather low to rather large values. As in [30], C rates of C = 
1/20, 1/10, 1/5, 1/2, 1,2,3,5,7, 10 were used, which ultimately aims at revealing 
rate capability of the underlying cell composition. 


The comparison of experiments and simulation results can be observed in 
Figure 4.8. Recall that the electrode structure with unmodified dense cathode 
particles is modeled by using the classical half-cell model and is denoted by 
p-NCM. Moreover, electrodes built-up by nanostructured secondary particles 
are modeled by means of the hierarchically structured half-cell model and 
are referred to as n-NCM-F850 and n-NCM-F900. The experimental values 
in Figure 4.8b were obtained by taking the mean value of the corresponding 
measurements in [30]. Notably, n-NCM-F850 and n-NCM-F900 show better rate 
capabilities as compared to p-NCM. Clearly, the same behavior can be observed 
in Figure 4.8a, where the simulation results also predict better performance for 
the nanostructured particle electrodes. Especially for rather high C rates of 7C 
and 10C, the hierarchically structured half-cell model results are very close to 


-l and 


the measured values. While the simulations estimate around 70mA hg 
50mAhg™! retained specific capacity for these C rates, experiments measure a 
little bit higher values of around 75 mA hg"! and 55mAhg~!. Note that in case 
of the classical half-cell model, the simulated value is a little bit higher in case of 
10C than the measurements. The simulation predicts approximately 30 mA hg !, 
while experiments found approximately 25 mAh g^!. Moreover, in case of both 
the classical and hierarchically structured half-cell model the retained specific 
capacity only decreases gradually for C rates between !/20C and 3C. Apparently, 
the simulation curves show a distinct drop between 3C and 7C, whereas the 
experimental curves decrease steadily from low C rates to higher C rates. 
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Figure 4.8: Rate capabilities of p-NCM, n-NCM-F850 and n-NCM-F900. a) Simulation results. b) 
Experiments from [30]. 


Given that the results clearly show a qualitatively good agreement with 
experiments, the proposed hierarchically structured half-cell model can be 
regarded as validated in a sense that it allows for investigations to predict at least 
qualitatively the influence of geometry, structure and material on electrochemical 
performance. Moreover, as is shown in the following, the model offers detailed 
information on how those parameters affect the cell quantities. For instance, 
the drop behavior mentioned before can be understood by investigating the 
local distribution of DOD, denoted here as dod. It is calculated similar to 
Equation (4.45), however, the volume-averaged lithium concentration of the 
solid phase at time ¢ is not evaluated for the whole cell, rather, the average 
concentration is computed at every point x along the cathode on cell level and 
at every point y along the secondary particle level. 


On cell level, the local distribution is computed as 


C,(f,X) — C. ius 
dod(t,x) = ER) "Cs (4.127) 
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The volume-averaged lithium concentration of the solid phase at time t at the 
position x is calculated as 


Тес 


T (t,x) = 3 y2c,(t,x) dy (4.128) 


Гес 0 


in case of the classical half-cell model, whereas this value is 


3 "sec 3 "prim Я 
ар у | ce Gandy, (4.129) 
Гес 0 0 


r prim 
in case of the hierarchically structured half-cell model. 


Figure 4.9 shows dod over the normalized cathode position for both the classical 
and hierarchically structured half-cell model. r is set as the end time of 
the simulation, which is determined by the cut-off voltage 3 V. The cathode 
positions (x — L*P)/LP9$ = 0 and (x — L**P)/LP^* = 1 refer to the separator- 
cathode interface and the cathode-current-collector interface, respectively, while 
(x — LP) /LP°S = 0.5 is the cathode middle position. Additionally, in Figure 4.9a 
- €, the state of charge distribution is shown for different C rates, namely 3C, 5C 
and 7C. 


It can be observed that in case of 3C, the distribution is constant over the 
cathode thickness. This applies to p-NMC as well as n-NMC-F850 and n-NMC- 
F900. While dod is almost 1 in case of n- NMC-F850/n-NMC-F900, the depth 
of discharge is lower in case of p-NMC. At 5C, all curves show decreasing 
behavior towards the cathode-current-collector interface. In case of n-NMC- 
F850 and n-NMC-F900, a small gradient is visible in the first half of the cathode 
whereas in the second half the gradient is more pronounced. In case of p- 
NMC, the decreasing behavior can be observed directly at the separator-cathode- 
interface. Overall, the dod values of the hierarchically structured half-cells are 
constantly larger than 0.60 whereas the classical half-cell model ranges around 
0.45. Finally, at 7C, the depth of discharge distribution of n-NMC-F850 and 
n-NMC-F900 looks similar to the one of p-NMC at 5C. However, for the 
nanostructured electrodes, a gradient is visible right at the separator-cathode 
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interface and the curve flatten towards the cathode-current collector interface. 
Moreover, the dod level of n-NMC-F850 and n-NMC-F900 is around the same 
value as well. The distribution in case of p-NMC is almost constant over the 
cathode region and is around 0.25. 


By comparing the dod distribution of the different cells, an important aspect 
considering cathodes with nanostructured secondary particles can be extracted. 
For increasing C rates, the cathode active material is more efficiently utilized 
in case of hieararchically structured cathodes, which yields to superior rate 
performance as compared to cathodes with monolithic secondary particles. In 
addition to the dod distribution along the x-direction, the depth of discharge 
distribution can be investigated along the radius in y-direction inside the 
secondary particle, which is referred to as dod" (7, y). 


On secondary particle level, the depth of discharge distribution along the radial 
position is calculated as 


Cre (t E y) m Cs init 


dod***(r, у) = 
od (t, y) = 


(4.130) 


Cs max — Cs init 


The volume-averaged lithium concentration of the solid phase at time 7 at the 
position y is calculated as 
hy); (4.131) 


in case of the classical half-cell model, whereas this value is 


3 


sec . 
cs (t,y) = — |, zd me y) dz. (4.132) 


4 prim 


in case of the hierarchically structured half-cell model. 
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Figure 4.9: Local depth of discharge along cell direction at different C rates. a) 3C. b) 5C. c) 7C. 
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In Figure 4.10, dod*** (7, у) is plotted versus both the normalized x and y position 
inside the secondary particle. The cathode position (x — LP) /LP95 = 0 refers 
to a secondary particle sitting at the separator-cathode interface, whereas (x — 
Г^®Р) /ІРӘ% = 0.5 and (x — L'*P)/LP^5 = | are secondary particles in the middle 
of the electrode and at the cathode-current-collector interface, respectively. 
Additionally, y/r,.. = 1 denotes to the secondary particle surface, whereas 
y/r4«, = 0 is the center. Again, the depth of discharge distribution is shown for 
the distinct C rates of 3C, 5C and 7C, respectively, and 7 is set as the simulation 
time at the cut-off voltage of 3 V. 


In general, the curves corresponding to n-NCM-F850 and n-NCM-F900 show 
less pronounced gradients as compared to the p-NCM curves. At ЗС, it can 
be observed that the surfaces of the secondary particles in both classical and 
hierarchically structured half-cell model show high states of discharge at around 
1.0. However, from the surface to the center of the secondary particles a 
more drastic decrease can be seen in case of p-NMC as compared to n-NMC- 
F850 and n-NMC-F900, where the latter two of these remain stable at a high 
level. At higher C rate of 5C, the differences in dod°““ distribution in the 
secondary particles become more visible. While n-NMC-F850 and n-NCM- 
F900 show a rather constant distribution, a clear gradient can be observed for 
p-NCM. Additionally, following the observations made for Figure 4.9, the state 
of discharge distributions at the secondary particle surfaces along the cathode 
direction decline from the separator-cathode interface to the cathode-current- 
collector interface. An interesting observation can be made for the case of 7C. 
Again, the gradients in the dod*** distributions inside the secondary particles are 
more pronounced for p-NCM than in the cases of n- NMC-F850 and n-NMC- 
F900. However, the depth of discharge at the secondary particle surface is 
actually higher at the separator-cathode interfaces in case of p-NMC than in the 
cases of n-NMC-F850 and n-NCM-F900. 


From the above observations of the dod**^ distribution, the conclusion from 
above can be extended. In addition to the enhanced utilization of the cathode's 
secondary particles due to an internal porosity, the depth of discharge distribution 
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Figure 4.10: Local depth of discharge along cell and secondary particle direction at different С rates. 


a) 3C. b) 5C. с) 7C. 
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and, thus, the concentration level of the active material is homogeneous. Even 
for larger C rates, the homogeneity remains stable inside the nanostructured 
secondary particles, which leads to a better rate performance as compared to 
the non-porous secondary particles. 


As was reported in [30], the investigated electrodes n-NNC-F850, n-NMC-F900 
and p-NMC show comparable loadings, i.e. active material mass per unit cross 
section. As was discussed above, the electrodes with the porous secondary 
particles show better utilization of the active material. This way, the energy 
density, i.e. energy per unit mass, is improved as well. 


From the above observations, the superior performance of hierarchically 
structured electrodes as compared to classical electrodes can be reproduced by 
means of the presented models. This means, for instance, the hierarchically 
structured half-cell model can be used to estimate sensitivity of geometry, 
structural and material parameters, which is the topic in Section 5.3. 


4.4 Summary 


In this chapter, a model for half-cells with hierarchically structured cathodes with 
porous secondary particles was proposed. As a first step, the well-established 
half-cell model based on the contributions of Newman and coworkers was 
recalled in detail. By means of the volume average method, the hierarchically 
structured half-cell model was proposed as a consistent extension of the classical 
half-cell model. The mathematical boundary conditions were acquired in detail 
and physically motivated for both models. Finally, the hierarchically structured 
half-cell model was qualitatively validated by comparing simulated and measured 
electrochemical performance based on real-world cathode structures taken from 
literature. The model successfully predicted the experimentally observed 
superior electrochemical performance of cathodes with porous secondary 
particles as compared to monolithic secondary particles. Additionally, the model 
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allowed for the investigation of this observation by discussing the local lithium 
concentration distribution of the solid phase of the electrode and secondary 
particle level, respectively. It was shown that nanostructured secondary particles 
differ from classical ones by having more of a homogeneous concentration 
distribution for higher C rates. This way, the available active material capacity 
could be better exploited leading to a higher performance. 
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In the following, the applicability of the in Sections 3 and 4 presented methods 
is demonstrated. On the one hand, the resistor network method is employed 
to derive prediction formulas for the effective resistance of porous secondary 
particles and for the effective transport properties of sphere packings. In both 
cases, transport is considered via the surface and the volume of the particles. 
On the other hand, the hierarchically structured half-cell model is used for 
parameter studies in order to reveal the influence of different structural and 
material parameters on the electrochemical performance of cathodes with porous 
active material particles. 


5.1 Effective transport properties of 
hierarchically structured spheres 


In Section 3.1, it was shown that the resistor network approach can be utilized to 
compute effective transport properties of sphere assemblies. As for classical cell 
models of the type presented in Section 4.1, for instance, this approach can be 
used to deliver effective transport properties of the porous electrode structure. 
That is to say, if the electrode is modeled as sphere packings. Concerning 
the hierarchically structured cell model from Section 4.2, the active material is 
modeled as spherical secondary particles built up by smaller spherical primary 
particles. In this case, in order to model effective transport properties through 
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the active material phase, certain parts of the resistor network method need to be 
adjusted accordingly. 


Consider the exemplary network of porous secondary particles, in Figure 5.1а. 
In Figure 5.1b, the equivalent node resistor network can be observed, where the 
secondary particle centers are the nodes and the contacts between the particles 
are identified as the resistors. In Section 3.1, it was pointed out that the crucial 
part in the framework of the RN is the computation of the resistances of the 
individual contact pairs. The most practical way is to have analytical formulas 
as in Equations (3.3) and (3.11), where the resistances can be evaluated based on 
the geometry of the overlapping particles. 


Therefore, the following section aims at deriving a formula describing the 
resistance between contact pairs of hierarchically structured spheres by their 
geometry and inner structure. Moreover, two cases will be considered. In the 
first case, the transport is through the volume of the primary particles and, in the 
second case, the transport is via the surfaces of the primary particles. 


Figure 5.1: Resistor network method for hierarchically structured spheres. a) Assembly of 
hierarchically structured spheres. b) Equivalent network with nodes at the centers of the 
spheres and resistances assigned to the contact pairs. 
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Model 


To begin with, similar to Equation (3.11), the resistance between two contacting 
hierarchically structured spheres is modeled as a series connection of two 
resistances. Moreover, those two resistances are computed by utilizing a 
model where each of the spheres are put between two conducting plates. In 
Figure 5.2, the replacement model of a hierarchically structured sphere between 
two conducting plates is sketched. The larger secondary particle comprises 
of smaller primary particles. The respective primary particle and secondary 
particle radii are ЕА 
particle is positioned between two conducting plates with a gap distance of 


and r,.., respectively. Furthermore, the secondary 


h**. which results in a contact radius of Fece. The internal contact radii 
Bur are between the primary particles with radii уйт апа P Clearly, 


the effective resistance of the model in Figure 5.2 is depends on all the above 
presented structural properties. A formula representing effective resistances of 
hierarchically structured spheres must therefore contain all those quantities. 


nsec 


Figure 5.2: Model of hierarchically structured sphere between two conducting plates. 


The approach chosen here was to conduct a variety of numerical experiments 
and deduce structural dependencies. The resistor network method was chosen to 
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perform the numerical experiments. However, before running parametric studies, 
the applicability of the RN on hierarchically structured secondary particles with 
inner porosity has to be verified. 


Verification of the RN for hierarchically 
structured spheres 


Similar to rectangular domains in Section 3.1, hierarchically structured spheres 
were generated, the structure was densified and the effective resistance was 
calculated using the resistor network method. Finally, the results were compared 
to finite element simulations. Concerning the generation of the assemblies, 
the RCP was modified to create randomly generated, overlap-free and densely 
packed collections of spheres inside a spherical domain. The initial structures 
were then further densified using the numerical sintering algorithm. In this 
case, the primary particles were fixed in space and their radii were successively 
increased until a certain densification criterion was reached. This way, the 
spherical shape of the secondary particle was maintained. 


Importing the data, setting up the models and evaluating the effective transport 
properties for both cases of volume and surface transport followed the 
descriptions from Sections 3.1.4 and 3.1.2, respectively. In order to model 
the conducting plates, in case of the RN, sphere caps of a certain height were 
removed from the secondary sphere, see Figure 5.3a, on two opposing sides in 
transport direction, see Figure 5.3b. Consequently, all primary particles, with 
their centers lying inside those virtual caps, were deleted, see the highlighted 
spheres in Figure 5.3b. Finally, the bottoms of the caps were chosen to be the 
conducting plates connected to the intersecting primary particles, see Figure 5.3c. 
This is where the boundary nodes were placed following the description in 
Section 3.1. 
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Figure 5.3: Hierarchically structured sphere between two conducting plates. 


In the finite element model, the primary particles were simply cut off at a defined 
position in the transport direction. This can be seen in Figure 5.4b and d. 
Concerning the FEM models, all nodes on one of the conducting plate surfaces 
were set to the same temperature, while the potential of the nodes on the opposite 
surface were dropped by AT — 1 with respect to this potential. Due to the 
mathematical equivalence of the transport problems considered in this work, the 
methods used for heat transport problems solved by finite element methods can 
be compared equivalently to electric transport problems solved by the resistor 
network method. Therefore, concerning the RN, an electric potential drop of 
AQ — 1 was imposed on the boundary nodes. Analogous to Equations (3.5) 
and (3.7), the effective resistance of a secondary particle was computed by 


AT А 
REM I and RRN = 2 (5.1) 
q eff 


respectively. In the FEM model, the resulting heat flux А EM was obtained at one 
of the surfaces with the applied temperature and, in case of the RN, the effective 
current Г was calculated using the resistor network solving scheme according 
to Section 2.3. 
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Figure 5.4: Verification models for the effective conductivity of hierarchically structured spheres 
using the resistor network and finite element method. Top: Transport via the surface, 
i.e. solid-surface. Bottom: Transport via the volume, i.e. solid-volume 


Test cases 


In the following, 15 scenarios of assemblies were created, according to the 
method described above. The scenarios were divided into three types of 
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assemblies, see Table 5.1. For all assembly types, the mean radius mean and 
the mean contact angle Omean were fixed to 1 and 15°, respectively. The only 
varying parameter was the standard deviation rg ranging from 0 to 0.25 such 
that, with increasing assembly type number, the polydispersity increased. The 
resulting volume fraction of the solid phase фона, теап and pore phase ópore mean 
as well as the mean contact radius rc,mean is presented in Table 5.1 as well. In 
addition, the transport is considered either through the volume or the surface of 
the primary particles. In the latter case, the surface is modeled by shells with a 
thickness of 0.1. r 


prim? 


i.e. 10% of the primary particle radii. 


Table 5.1: Structural parameters of secondary particles for the verification of the transport through 
the volume and via the surface of solid spheres. 


Type mean ro Omean Ösolid,mean Öpore,mean l'c,mean 

1.00 0.00 15 0.653 0.327 0.257 
2 1.00 0.10 15 0.650 0.350 0.245 
3 1.00 0.20 15 0.646 0.354 0.234 


Evaluation and discussion 


The results of this investigation can be observed in Figure 5.5. Effective 
resistances REN and REEM are calculated by Equations (3.8) and (3.17), 
respectively, and normalized by the bulk or shell conductivity. On the horizontal 
axis the results provided by the FEM analysis are shown and on the vertical axis 
the corresponding values by the RN are drawn. The black solid line in this figure 
represents a perfect match of the two results whereas values above or below 
indicate an over- or underestimation of the effective transport property by the 
RN, respectively. 
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Figure 5.5: Verification of the resistor network method for hierarchically structured spheres using the 
finite element method. a) Transport through the volume of primary particles. b) Transport 
via the surface of primary particles. 


It can be seen that, while the mean error is around @ = 6% in the case of 
surface transport, the mean error is almost zero in the case of volume transport. 
Therefore, given the considerably low errors, the above observations allow the 
conclusion that the resistor network approach is suitable to compute the effective 
transport properties of hierarchically structured spheres between two conducting 
plates. 
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Parametric studies and prediction formulas 


In order to model the effective resistance, the prediction formula 


Porimym : sec ‚ primym | 1 


ofi 
Ке = Ё". * Poulk,shell (5.2) 


Porim,cym Porimyn Sm sec 


is proposed, where At is a fit function and the rest of the parameters are 
interpreted as follows. Here, the first model parameter is the ratio ”primm/r 


prim,c,m? 
which is interpreted as the degree of densification of the secondary particle and 


relates the mean primary particle radius r to the mean contact radius of 


prim,m 


the primary particles r Furthermore, 7/7... describes the tortuosity 


prim,c,m* 
effect and takes the secondary particle radius г into account. Finally, "sec,c/rsec 


is the ratio of the contact radius of the conducting plates r, to the secondary 


sec,c 
particle radius and рь she is either the bulk or the shell, resistivity, depending 


on if volume or surface transport is considered. Additionally, sm/r relates the 


prim,m 


shell thickness Sm to the mean primary particle radius, where, in case of volume 
transport, this parameter is omitted. 


The degree of densification "prim,m/r can theoretically range within the limits 


rim,c,m 
of 1 and ce. In the first case, the dontct radii are of the same size as the primary 
particles, which means that the structure can be regarded as bulk material. The 
second case represents a structure where the primary particles have no contact 
to each other, which results in no conducting pathways. The tortuosity effect 
described by "sec/r 


closer the primary and secondary particle sizes are, the shorter the conducting 


primm 18 Within the limits of 1 and ©. Roughly speaking, the 
paths within the secondary particle. As a consequence, the effective resistance 
decreases. On the other hand, the effective resistance increases if comparably 
small primary particles form more tortuous and longer pathways. Finally, the 
values of ”primm/sm are within the limits of 1 and eo. The latter of which resembles 
a mean shell thickness of 0 and no transport paths via the surface. Consequently, 
the effective resistance tends to infinity as well. In theory, a shell thickness of 
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the size of s,, provides the best transport paths via the surface and, therefore, the 
lowest effective resistance. 


In the following, the verification of the prediction formula in Equation (5.2) 
is demonstrated by 100 randomly generating hierarchical sphere structures 
following the tool chain of RCP (modified to account for spherical domains), 
numerical sintering and RN, see Section 3.1. For each of the scenarios, the 
combination of the model parameters were chosen at random within the limits 
shown in Table 5.3. 


Table 5.3: Variation of parameters of the hierarchically structured spheres. 


Parameter Value 

prim m T sim c m 2.4 +. 1 1.5 
Tsec] rim m 13.4...41.0 
ta 0.01 n 0.09 
ест 0.2. . .0.4 
Pbulk,shell 0.001...0.1 


Figure 5.6 shows the results of the numerical experiments. In Figure 5.6a, 
transport through the volume and, in Figure 5.6b, transport via the surface is 
shown. The effective resistance is plotted over the dimensionless ratio "sec.c/r im m’ 
It can be observed that, on the one hand, the results can be represented in 
a dimensionless form using the parameters in Table 5.3. On the other hand, 
the curves approximately lie on a simple rational !/x-type function, which only 


depends on the ratio "sec,c/r, This function can be given by 


prim,m* 


ee (5.3) 
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where ag; 25 6, in case of volume transport, and ag œ~ 20, in case of surface 
transport. 


Using the prediction formula from Equation (5.2), the behavior of the effective 
resistance in terms of the structural properties can be investigated. For instance, 
in both the volume and the surface transport case, the effective resistance tends 
to increase to rather large numbers if "sec,c/r 


primm 


approaches small values. This 


means that, given a constant г. ., for increasing r 


sume the number of primary 


particles in contact with the conducting plate decreases leading to higher effective 


resistances. Vice versa, an increasing ratio Of ‘see,c/r leads to a smaller 


prim,m 


effective resistance, which means that more primary particles are in contact with 
the conducting plates. 


pRN,surf 
Кот 


Riitsurf 
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Figure 5.6: Prediction formulas to estimate the effective resistance of hierarchically structured 
spheres between two conducting plates. a) Transport is via the volume of the primary 
particles. b) Transport is via the surface of the primary particles. 
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Conclusion 


Equation (5.2) represents the effective resistance of a hierarchically structured 
secondary sphere between two plates. The secondary sphere is composed of 
smaller primary spheres. The formula can be used to model transport either 
through the volume or via the surface of the primary spheres. Moreover, it 
can be used to model transport between two contacting hierarchically structured 
secondary particles. Practically, and similar to Equations (3.3) and (3.11), the fit 
formula can be used in the framework of the resistor network method to calculate 
the effective conductivity of an assembly of hierarchically structured spheres. 


5.2 Effective volume and surface transport 
properties of sphere assemblies 


The knowledge of effective transport properties in granular materials such 
as battery electrodes [18, 150] plays an important role for battery modeling. 
Usually, effective transport properties of granular materials are approximated 
by the volume fraction of the transport phase under consideration, see [69] for 
a review. Battery modelers tend to use a Bruggeman type relation [71, 138] 
to estimate effective transport properties of the porous electrodes [36, 41, 46]. 
In [92], the resistor network method was used to show that this type of empirical 
relation works well for the pore phase of sphere assemblies with a polydisperse 
size-distribution. In order to model effective transport properties of particle 
networks, however, the Bruggeman relation is not well suited [67]. 


In order to overcome this drawback, adequate prediction formulas need to 
be found. Therefore, in the following, the practical use of the solid-volume 
model from Section 3.1.2 and the solid-surface model from Section 3.1.4 is 
demonstrated. The models are used to develop prediction formulas to compute 
effective transport properties of assemblies of overlapping and equal-sized 
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spheres, where the transport can be through the volume or the surface of the 
spheres. 


Model parametrization 


A series of 100 scenarios of randomly distributed and monosized sphere 
assemblies were created in a box-shaped simulation domain using a certain set 
of model parameters. The model parameters were chosen as follows. A porosity 
range pore between 0.1 and 0.4 was targeted representing, on the one hand, a 
rather dense packing and, on the other hand, one which is near the percolation 
threshold of monosized spheres [76]. Moreover, the sphere radii r, were between 
0.1 and 0.8 representing an arbitrary but rather large range of particle sizes. This 
particle size range is comparable to the primary particle sizes possible in [30], 
which ranged between 0.179 and 0.6um. Note, however, an arbitrary length unit 
lu is chosen, since the investigated transport parameters are dimensionless and 
independent of the chosen length unit. Finally, in case of surface transport, a 
constant shell thickness was set to 0.01. 


Analogous to Section 3.1.2, the initial sphere packings were generated using 
the RCP. The structures were then further densified by means of the numerical 
sintering algorithm. Additionally, the numerical sintering algorithm was slightly 
modified. In contrast to the algorithm described in Section 3.1.2, the box- 
shaped simulation domain was isotropically and successively reduced such that 
the spheres inside the domain were moved upon each other. This was done until 
a densification criterion, i.e. the targeted porosity, was achieved. In this way, the 
particle sizes could be preserved. 


157 


5 Applications 


Results and discussion 


The effective transport parameters through the volume and via the surface of 
the spheres were calculated using the resistor network method as presented in 
Section 3.1.2 and Section 3.1.4, respectively, and evaluated in Figure 5.7. Here 
the effective transport parameters are drawn over the porosities and particle radii. 
In Figure 5.7a, the effective transport parameters of the volume transport case 
is shown whereas, in Figure 5.7b, the surface transport case can be observed. 
Recall from the above mentioned sections that the effective transport parameters 
are defined as the effective transport properties normalized by the transport 
coefficients of the particle volumes or shells. In the present investigation, without 
loss of generality, they both set to unity. 


Apparently, in Figure 5.7a, the effective volume transport parameters only 
depend on the porosity of the system. For high porosity values the effective 
transport parameters are low while for low porosities they increase. Therefore, 
the numerical results are fitted by the bi-linear surface fit 


kf = а) +a, pore +aprp, (5.4) 


where the parameter values are ag = 1.0, ај = —2.130 and a» = 0.01. Here, ag 
is chosen such that, for the special case of a completely dense structure, i.e. if 
Ópore — О and rp — 0, the effective transport parameter becomes 1 resembling 
bulk material properties. 


In Figure 5.7b, the effective transport parameter via the surface is shown. 
Remarkably, the effective transport parameter is not only dependent on the 
porosity but also on the particle radius. The best effective transport parameter 
is achieved for low particle radii and low porosities. Additionally, it can be 
observed that for a constant porosity, the effective transport parameter can be 
enhanced by decreasing the particle radius. Finally, the resulting values where 


fitted by 
jfit 1 1 


Коке . ; 
id (1 + pore)” (14 rg)^i 


(5.5) 
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where the fitting parameter values are bọ = 6.717 and b, = 5.332. Note that 
the function was chosen such that, for the completely dense structure case, i.e. 
Ópore — О and ry > 0, Equation (5.5) becomes bulk material property. 
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Figure 5.7: Evaluation of transport parameters by the parametric study. a) Transport through volume. 
b) Transport via surface. c) Comparison of surface and volume transport over the porosity. 


Consider the special case of linearly increasing particle radius and decreasing 
porosity. A similar case was observed in [151]. There, pellets were 
investigated made of conductive particles which were pressed together and 
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sintered afterwards. By adjusting the pressure and sinter temperature, different 
particle radii and the pellet porosities were achieved leading to increasing particle 


radii while, at the same time, decreasing porosity. In Figures 5.7a and 5.7b, such 
a path is indicated by pus or 
respectively. In Figure 5.7c, Equations (5.4) and (5.5) are evaluated along this 


path. In particular, the curves are plotted over porosity. While LU m 


in case of volume or surface transport, 


decreases 
with increasing porosity in the volume transport case, which can be expected, 


7fit,path 7fit,path 


the surface transport curve k „p гіѕеѕ. Clearly, k į decreases because of 


reducing overlaps of particles. However, reducing overlaps apply also to pae 
Therefore, this effect seems to be compensated by decreasing particle sizes 


because, due to that, the specific surface area increases as well. 


Apparently, the specific surface area curve along the given path is similar to 
the evolution of the transport parameter via the surface. Given that the specific 
surface area is intrinsically accounted for by the transport angle of the resistance 
in Equation (3.11), this directly effects the overall effective conductivity. 


Total free 
surface area 


Free surface area 


Sphere cap 


Figure 5.8: Computation of specific surface area for assemblies of spheres. 


The specific surface area is defined as the total free surface area divided by the 
simulation box volume. The computation of the specific surface area—as used in 
this work—is sketched in Figure 5.8. For a given assembly of spheres depicted 
on the left-hand side, the total free surface area is calculated by summing up the 
surface areas of each individual particle. The magnification on the right-hand 
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side exemplarily highlights a single particle out of the assembly. The individual 
free surface area is computed by subtracting the surface area of the sphere caps, 
which are formed by the overlaps with the neighboring particles, from the total 


surface area of the particle. 


In Figure 5.9, the evolution of specific surface area for the present investigation 


can be observed. In Figure 5.9a, the resulting values are fitted by 
fit co c2 


AA (1 + pore)?! | (1 +ту)% u (5.6) 


where the fitting parameters аге со = 3.636, cy = —1.104, c2 = 5.196 and 
сз = 5.397. The evolution of specific surface area along the given path Qm 


can be seen in Figure 5.9b. 
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Figure 5.9: Specific surface area by the parametric study. 


161 


5 Applications 


Conclusion 


The investigation results reveal that Equations (5.4) and (5.5) can be used to 
model the effective transport parameters via volume and surface of mono-sized 
sphere assemblies with varying porosity and particle sizes. In case of cell 
modeling, Equation (5.4) can be used to calculate the transport properties of 
the active material phase, given that it is modeled by spheres. Additionally, 
since in Equation (5.5) thin layers of shells are assumed, it can be utilized to 
model effective transport properties of high-conductive carbon coatings of active 
material. Finally, the above findings support the observations in [151]. There, an 
increase in effective ionic conductivity was observed in case of larger porosities. 
This increase was attributed to the dominant transport via the surfaces of the 
particles due to an increase in surface area. 


5.3 Parameter study on 
hierarchically structured 
electrode performance 


In Section 4.3, the hierarchically structured half-cell model, as proposed in this 
work, was validated qualitatively by experiments taken from literature. It was 
concluded that it can be used to investigate the influence of the electrode’s 
geometry, structure and material on the cell’s electrochemical performance. 
By means of this model, the difference between classical and hierarchically 
structured cathodes in terms of sensitivity to the diffusion coefficient and 
electronic conductivity of the active material was investigated in [137]. In the 
following study, however, the focus is on the hierarchically structured cathodes 
only and a larger set of parameters is taken into account. The influence of 
diffusion coefficient and electronic conductivity of the active material as well as 
secondary particle morphology on the rate capability is analyzed. To this end, the 
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cathode structure n-NMC-F900 from Section 4.3 was used as a reference and the 
mentioned parameters were varied as follows. Note that, if not stated otherwise, 
the other parameters concerning geometry, structure, transport and material were 
taken from Tables 4.1, 4.2, 4.3 and 4.4. 


The electronic conductivity was assumed between 1-1073, 1-1074 and 
1-107? S m^! representing a rather large spectrum of high and low values 
and corresponding to lithiation states of y — 0.9, 0.95 and 1.0 in 
Li4Ni;/3Mn,/3Co1/3O»?, see Equation (4.124). The diffusion coefficient was set 
as 1-107, 1-107!4 and 1-107 PD m?s-, which are the extreme values found 
in literature [142-146]. The secondary particle morphology is characterized by 
the porosity, primary particle radii, specific surface area and effective transport 
parameter. In this investigation, the primary particle radii and effective transport 
parameters were varied while keeping the porosity fixed. Note that, indirectly, the 
specific surface area is influenced by changing the radii since it is calculated as 
а = 3e [136]. In order to study the influence of the transport via the primary 
particle network, the effective transport parameter 42% is varied between 
values of 0.25, 0.50 and 0.75, where the first and the middle value is comparable 
to the lowest and largest parameters which were found for n- NMC-F850 and 
n-NMC-F900, respectively. An additional effective transport parameter of 0.75 
is chosen to see if an improvement of transport path quality would influence 
performance. 


For all 81 combinations of the above described parameters, the rate performance 
was simulated by using C rates from 1/20 to 10 and the simulation was stopped at 
a cut-off voltage of 3 V. Results of this investigation can be found in Figure 5.10, 
where some representative cases are selected to be shown. From the top to the 
bottom row, the electronic conductivity is 1: 1073, 1-1074 and 1-10? Sm-!, 
respectively. From the left to the right column, the transport coefficient, the 
primary particle radii and the diffusion coefficient is varied. 


In general, it can be observed that from the top to the bottom the rate capability 
decreases for decreasing electronic conductivity. While a distinct specific 
capacity drop in case of the largest two conductivity values can be observed 
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after 3C, the performance drops after 1C for the case of lowest conductivity. 
Note, however, that the gradient of the drop is steeper in case of the larger two 
conductivity values. On the other hand, for the second and third column, the 
variation of particle radii and diffusion coefficient does not seem to influence 
the shape of the performance curve as the curves within one diagram are almost 
identical. However, in case of the first column, where the transport coefficient 
was varied, a clear difference between the curves can be observed for decreasing 
conductivity of the active material. Especially, in the lower left diagram, it can 
Teon,sec 


be seen that for increasing values of Ке”, the rate performance decreases more 
visibly. 


From the above observations, some conclusions can be drawn. Varying the 
diffusion coefficient within the range presented here does not seem to influence 
the performance. While it was found in [137] that for classical cathodes the 
rate-limiting factor is the diffusion of lithium into the monolithic active material 
particles, in case of the hierarchically structured cathodes, sensitivity to the 
diffusion coefficient could not be observed, see also [53]. In addition to that, 
changing the primary particle radii, in the range presented here, has no visible 
effect on the performance. This indicates that the primary particle radii in 
combination with the diffusion coefficient still provide diffusion paths which 
are short enough to not have an impact. However, electronic conductivity of 
the active material has a strong effect on rate capability. On the one hand, this 
can be explained by assuming that lithium diffusion via the electrolyte phase 
of the secondary particle pores is favored. Here, the diffusion coefficients are 
orders of magnitudes higher than in the solid phase [146, 152]. On the other 
hand, low electronic conductivity values lead to a lower performance, which 
was also observed in [53, 152]. Therefore, electronic conductivity becomes rate- 
limiting. In addition to this observation, this investigation revealed the secondary 
particle morphology plays a crucial role when the electronic conductivity reaches 
lower values. In other words, poor inter-primary particle contact becomes more 
problematic in case of low electronic conductivity, which was also observed 
in [152]. 
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Figure 5.10: Parameter study rate performance. Varying electronic conductivity of active material 
from top to bottom row. Varying transport coefficient, primary particle radii and 
diffusion coefficient from left to right column. 


The effect of the transport coefficient on the electrochemical performance is 
further investigated with the help of Figure 5.11, where, analogous to Section 4.3, 
the depth of discharge distribution along the secondary particle radii is shown. 
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The lowest electronic conductivity of 1 - 10755 т! is chosen at a C rate of 3C 
as, in this case, the differences in performance are visible the most. Three cases 
of transport coefficients are discussed, where it can be observed that for all cases 
the value of dod°““ is around 1 at the surfaces of the secondary particles. From 
the surface to the particle center, the depth of discharge decreases slowly at first 
but at a certain point shows a distinct drop and decreases slowly again until the 
end. In case of higher values of £275, this distinct drop occurs more towards 
the surface of the secondary particles as well as results in lower dod°° towards 
the center. This leads to a lower depth of discharge DOD for the whole cell 


represented by lower retained specific capacity. 
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Figure 5.11: Depth of discharge distribution along secondary particle radius direction at 3C. 


Poor inter-particle contact and low electronic conductivity of the active material 
can explain the observations in Figure 5.12, which is based on investigations 
in [137] and where the two cathode structures n-NMC-F850 and n-NCM-F900 
are compared to each other. Experiments in Figure 4.8b show that for higher 
C rates, the two rate performance curves diverge from each other. Notably, 
however, n-NCM-F900 performs a little better than n-NCM-F850 even though 
the latter of which shows more preferable structural properties in terms of 
higher specific surface area of primary particles for electrochemical reactions 
inside the secondary particles and smaller primary as well as secondary particles 
leading to shorter diffusion paths. For the simulation, the rate performance for 
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the two structures were calculated using the geometry, structure and material 
parameters from Tables 4.1, 4.2, 4.3 and 4.4. However, the secondary particles’ 
electronic conductivity was set to the lowest value of 1-10°°Sm-! in both 
cases. This low value could be justified as a first simple measure to account 
for additional resistance effects such as contact resistances between primary 
particles beyond the purely geometrical bottleneck effects considered so far. The 
results show a clear difference in performance, where n-NMC-F900 achieved 
better rate capability than n-NMC-F850. By comparing the effective transport 
parameters from Table 4.3 of the solid phase for both cases, it can be seen 
that the value of {202% = 0.46 in case of n-NMC-F900 is more than twice 
as large as it is in case of n-NMC-F850, which is 0.20. As was pointed out 
before, effective transport parameters reflect the quality of transport paths which, 
ultimately, characterizes the connectivity of primary particles. Adding to the 
fact that electronic transport is predominantly via the primary particle network 
of the secondary particles [153], the combination of low electronic conductivity 
and poor connectivity of primary particles may lead to low electrochemical 
performance. Therefore, in view of the above observations, the rate-limiting 
factor is the effective electronic conductivity of the primary particle network 
inside the porous secondary particles. 
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Figure 5.12: Simulation results based on real electrode structure and using low electronic 
conductivity. 
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Conclusion 


A large scale parameter study was conducted by means of the qualitatively 
validated hierarchically structured half-cell model from Section 4.3. The 
influence of electronic conductivity, transport and diffusion coefficients of the 
active material as well as the primary particle radii on the rate performance was 
discussed. It was found that diffusion coefficients and primary particle sizes 
have barely visible effect as the diffusion paths might still be short enough 
to not have an impact on performance. However, it could be shown that low 
electronic conductivity and large transport coefficient influence the rate capability 
more significantly. It was concluded that the combination of low electronic 
conductivity and poor inter-connectivity of the primary particle network inside 
the secondary particles might become rate-limiting in case of hierarchically 
structured cathodes. 
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In the present work, the resistor network method was extended with regard to the 
transport through the solid and the pore phase of granular systems represented 
by sphere packings with polydisperse size-distributions. As for the solid phase, 
transport through the volume, via the surface or a mix of both were considered. 
For all cases, appropriate analytically derived formulas from literature describing 
resistance between two single particles were used or combined accordingly. 
Finally, these single contact models were embedded into the framework of the 
resistor network method in order to compute effective transport properties. All 
the proposed models—-single contact as well es effective transport models — were 
verified using finite element methods. 


Regarding the pore phase, a novel method concerning the computation of 
transport through the pore phase was developed. The key was the so-called 
Laguerre tessellation, where the pore phase of the system is decomposed into 
cells with each of them containing a single particle. Consequently, the cell 
nodes and edges were the basis of equivalent resistor networks. The nodes were 
identified as pore centers and the edges were the pore throats. As an extension, 
this model was modified to account for more than one conducting species in the 
pore phase. Both methods were either verified using the finite element method or 
validated with the help of experiments taken from literature. 


In the application part of this work, it was shown how the resistor network method 
for volume and surface transport can be used to conduct large scale parameter 
studies. Due to the resource and time efficiency of this method, a large variety of 
different structural compositions can be simulated in a short amount of time. In 
this work, the resulting database allowed for the derivation of prediction formulas 
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for the effective resistance of porous secondary particles and porous electrodes 
consisting of mono-sized and overlapping spheres. 


A model for half-cells with hierarchically structured cathodes with porous 
secondary particles was proposed. As a preparation step, the well-established 
half-cell model in accordance to the contributions of Newman and coworkers 
was recalled in detail first. It was pointed out that the mathematical framework 
of this model can be derived equivalently by employing volume averaging 
methods. By means of the same methods, the hierarchically structured half- 
cell model was proposed as a consistent extension of the classical half-cell 
model. For both models, the mathematical boundary conditions were acquired 
in detail and they were physically motivated. The hierarchically structured half- 
cell model was applied to real-world cathode structures taken from literature. 
In doing so, the model was validated qualitatively by comparing simulated 
and measured electrochemical performance as it successfully predicted the 
experimentally observed superior electrochemical performance of cathodes with 
porous secondary particles as compared to monolithic secondary particles. 
Additionally, it allowed for the investigation of this observation by discussing 
the local lithium concentration distribution of the solid phase along electrode 
and secondary particle direction. It was shown that nano-structured secondary 
particles differ from classical ones by having more of a homogeneous 
concentration distribution for higher C rates. This way, the available active 
material capacity could be better exploited leading to a higher performance. 
Finally, the hierarchically structured half-cell model was applied for parameter 
studies. By varying the diffusion coefficient and electronic conductivity of the 
active material as well as the secondary particle morphology, it was found that 
the combination of low conductivity and poor inter-connectivity of the primary 
particle network inside the secondary particles might become rate-limiting. 


As a possible future work, the resource and time efficiency of the resistor network 
as well as the presented half-cell models can be used for further investigations. A 
large variety of artificial particle-based structures can be generated and effective 
transport properties for the solid and pore phase can be computed representing 
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the solid and the electrolyte phase in the cell models. These parameters can 
be incorporated into the cell models in order to investigate the impact on the 
electrochemical performance. This way, the electrochemical performance for a 
large set of parameters regarding the geometry, structure and material properties 
can be computed efficiently and, ultimately, the best combination can be chosen 
to design optimized electrodes. 
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List of symbols 


Constants 
Name Description Value Unit 
Faraday’s 1 
F 96485 Cmol 
constant 
R DER 8.314 Jmol-! K-! 
constant 
A d 
Na “лз 6.02214076- 1023 то! 
constant 
Physics 
Name Description Unit 
Фф Electric potential у 
1 Electric current A 
i Electric current density Am? 
K Electric conductivity Sm! 
p Electric resistivity Qm 
Psurf Electric surface resistivity Q 
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List of symbols 


Physics cont. 


Name Description Unit 
Koon Electronic conductivity Sm! 
кїп Ionic conductivity Sm! 
Rx Electric resistance Q 
Gy Electric conductance S 

T Temperature K 
da Heat flux W 
q Heat flux density Wm^? 
À Thermal conductivity Wm K^ 
R, Thermal resistance KW! 
G, Thermal conductance WK"! 
c Mass concentration mol m? 
D Mass diffusion coefficient m?^s^! 
J Molar flux mols-! 
j Molar flux density molm~*s~! 
t Time S 
Chemistry 
Name Description Unit 
ft Mean molar activity coefficient [-] 
Signed stoichiometric coefficients of 
S4, S—, S0 [-] 
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cation, anion and solvent 


List of symbols 


Chemistry cont. 


V4, V 


ko 


io 


Description 


Unsigned stoichiometric coefficient of 
the cation and anion 


Signed charge numbers of cation and 
anion 


Transference numbers of cation and 
anion with respect to the solvent 


Concentration of the solvent 
Number of moles of a species 


Asymmetric charge-transfer 
coefficient 0 < а < 1 


Effective reaction-rate constant mol? ^ 


Exchange current density 


Resistor network method 


Unit 


m^-3«g-1 


Name Description 

F Flux 

R Resistance 

p Potential 

"n Effective current 

Kost Effective conductivity 
Kpulk Bulk conductivity 
Knell Shell conductivity 
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List of symbols 


Resistor network method cont. 


p surf 
Ppore 
Laomain 


Adomain 
R 


R 
R 


solid, vol 
solid,surf 
solid,mix 
pore,vol 


pore,mix 
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Description 


Core conductivity 

Effective transport parameter 
Volume fraction of the pore phase 
Volume fraction of the solid phase 
Radius 

Contact radius 

Standard deviation of radii 
Contact angle 

Transport angle 

Shell thickness 

Surface resistivity 

Pore resistivity 

Domain length 

Domain cross section 
Solid-volume resistance 
Solid-surface resistance 
Solid-mix resistance 

Pore-volume resistance 


Pore-mix resistance 


List of symbols 


Cell model 


Name 


Znsec 


Ф», Po Ф; 


™|SEC 


Фе, Pe» 9, 


Description 


Electronic potential in the solid phase, 
on cell and on secondary particle level 


Ionic potential in the electrolyte 
phase, on cell and on secondary 
particle level 


Electric overpotential, on cell and on 
secondary particle level 


Volume-averaged production term 
of electronic charge from inside the 
secondary particle level 


Volume-averaged production term of 
ionic charge from inside the secondary 
particle level 


Volume-averaged production term of 
lithium ions from inside the secondary 
particle level 


Concentration in the solid phase, on 
cell and on secondary particle level 


Concentration in the electrolyte phase, 
on cell and on secondary particle level 


Electronic conductivity in the solid 
phase 


Ionic conductivity in the electrolyte 
phase 


Diffusional conductivity in the 
electrolyte phase 


Unit 
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List of symbols 


Cell model cont. 


Name 


eon,sec 


eon 
Ky eft K, eff 


eon,pos 
seff 


ion ion,sec 
Ke etf K, eff 


ion,sep _ion,pos 
K, etf ? K eff 


ion ion,sec 
K 


Kp ett Kp eff 


ion,sep _ion,pos 
Kp eff > Kp eff 


sec 
D, ap Рет 
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Description 


Effective electronic conductivity in the 
solid phase on cell and on secondary 
particle level 


Effective electronic conductivity in the 
solid phase on cell level in the positive 
electrode region 


Effective ionic conductivity in the 
electrolyte phase on cell and on 
secondary particle level 


Effective ionic conductivity in the 
electrolyte phase on cell level in the 
separator and positive electrode region 


Effective diffusional conductivity in 
the electrolyte phase on cell and on 
secondary particle level 


Effective diffusional conductivity in 
the electrolyte phase on cell level in 
the separator and positive electrode 
region 

Diffusion coefficient in the solid phase 


Diffusion coefficient in the electrolyte 
phase 


Effective diffusion coefficient in 


the electrolyte phase on cell and on 
secondary particle level 


Unit 


List of symbols 


Cell model cont. 


Name 


sep pos 
D, еф De ett 


sec 
dse, Ase 


Se; sec 


J se 2, Ј ѕе 

ф ѕес 
S? Ts 

Фф sec 
е° re 


фу 


Q 
Be 


U, cell 


5 tot 
LSP, [pos [to 


Acell 
үроз 


l'sec» Porim 


Description 


Effective diffusion coefficient in the 
electrolyte phase on cell level in the 
separator and positive electrode region 


Specific surface area between solid 
and electrolyte phase on cell and on 
secondary particle 


Butler-Volmer type mass flux density 
at the solid and electrolyte interface 
on cell and on secondary particle 


Volume fraction of the solid phase on 
cell and on secondary particle level 


Volume fraction of the electrolyte 
phase on cell and on secondary 
particle level 


Volume fraction of the filler phase on 
cell level 


Ampere-hour capacity 
Equilibrium potential 
Cell voltage 


Separator, positive electrode and total 
thickness 


Cross section area of the cell 
Volume of the positive electrode 


Secondary particle and primary 
particle radius 


Unit 
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List of symbols 


Cell model cont. 


Name Description 


MAM Mass of the active material 


Gravimetric density of the active 


PAM . 
material 


Sub- and superscripts 


Кет 


Name Description 
init Initial value 
max Maximum value 
surf Surface 
eff Effective parameter 
ion Ionic 
eon Electronic 
$ Solid phase 
e Electrolyte phase 
pos Positive electrode 
sep Separator 
tot Total 
cell Cell 
sec Secondary particle 
prim Primary particle 
spec Specific 
avg Average 
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List of symbols 


Sub- and superscripts cont. 


Name Description 
res Resulting 
Acronyms 
Name Description 
AM Active material 
CB Carbon black 
RVE Representative volume element 
FEM Finite element method 
XFEM Extended finite element method 
RN Resistor network method 
LIB Lithium-ion battery 
DEM Discrete element method 
RCP Random close packing algorithm 
PDE Partial differential equation 
ODE Ordinary differential equation 
BC Boundary condition 
soc Local state of charge 
dod Local depth of discharge 
SOC Global state of charge 
DOD Global depth of discharge 
OCV Open circuit voltage 
NMC Nickel Manganese Cobald Oxide LiNi; /3Mg1/3C01/302 
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List of symbols 


Acronyms cont. 


Name 


EMT 
FIB/SEM 
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Description 


Effective medium theory 


Focused Ion Beam Scanning Electron Microscopy 


A Mathematical operations 


A.1 Divergence operator 


The divergence operator in the cartesian coordinate system is 


v.p 0% 95 ӘР, 


Ox ду д: ' кы 


where x,y,z are the cartesian coordinates and F,,F,,F, are the respective 
components of the flux vector. 


The divergence operator in the spherical coordinate system is 


- 101 1 д 1 д 
УВ ај — ыу ге 
P ( к) а rsin(®) dO (Bo) rsin(®) 99/9 к 


where r,O,d are the spherical coordinates and F,,Fe,Fe are the respective 
components of the flux vector. 
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A Mathematical operations 


A.2 Integration 


Cartesian domains 


The volume integral of a function f(x,y,z) is 


111 tæra = Г [тода (А.3) 


where L,,L,,L; is the edge length of a box. In case of the function being 
independent of y, and z, the volume integral is 


L pD op L 
/ / | f(x) dxdydz = DL? / f(x) dx (A.4) 
0 Jo Jo 0 

and the volume average is 


DL ft fade _ 1 (P 


= dx . (A.5) 
V pj € 
Spherical domains 
The volume integral of a function (л Ө,Ф) is 
2л qx pR 
f f | f(r,0,®)r’sin(®)drd@d® . (A.6) 
0 Jo Jo 


In case of spherically symmetric conditions, the angular dependencies are 
constant such that the concerning volume integrals reduce to 


An / y fr)rdr . (A.7) 
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A.2 Integration 


Furthermore, the volume average of the above equation is 


An[2of(r)hdr 3 ү 
DT = =. > (ғ) ан. (А.8) 
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B Physics 


B.1 Specific capacity from stoichiometry 


By using Faraday’s law, which states that the rate of production of a species is 
proportional to the current, and the total mass produced is proportional to the 
amount of charge passed multiplied by the equivalent weight of the species [47]: 


> s;M;lt 


m B.1 
oe nF ’ 8:0 


where m; is the mass of species і produced by a reaction in which its 
stoichiometric coefficient is s; and n electrons are transferred, M; is the molar 
mass, F is Faraday’s constant, and the total amount of charge passed is equal to 
the current / multiplied by time т. Rearranging Equation (B.1) brings 


Tt nF 

mi siMi (B.2) 
— 

Ог рес 


It 
т? 


which yields the specific capacity of species i as О; spec = where It = Qi is 


taken as the total capacity of species i. 


The above mentioned coefficients of species i refer to the electrode reaction 


s_M_* +84 Му2+ cs) Мо — ne”, (B.3) 
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B Physics 


where M is the symbol for the chemical formula of species i, and z_, 2; are 
signed charge numbers of cation and anion [10]. 


As for the positive electrode, the reaction in a lithium-ion cell is 
M+Lit +e === ШЫМ, (B.4) 


where the M stands for some metal oxide in the positive electrode material [10]. 
Rearranging Equation (B.4) into the form of Equation (B.3) as 


LiM-M-Li* =e", (B.5) 


leads to the coefficients z іф = 1, 5 0,5, = —1, so = О апаи = 1. 


In the following, Equation (В.б) is used to calculate the specific capacity of 
the positive electrode active material LiNi;/3Mnj/3Co 1/302, which is used in 
this work. Therefore, the molar mass М; = Mini, /3Mnı /3C01/302 = Mwwc is 
computed as the stoichiometrically weighted sums of each element using their 
atomic weights. From Table В.І, the sum molar mass is 96.46 g mol" !. 


Table B.1: Atomic weights of the elements in Li Ni; /3 Мп /3C01 /302. 


Dl. Atomic weight Total weight 
Element Stoichiometry 21 Ži 
[g mol] [gmol”] 

Li 1 6.941 6.941 

Ni 13 58.693 19.564 
Mn 13 54.938 18.313 
Со 13 58.933 19.644 
О 2 15.999 31.998 


208 


B.2 Density from X-ray data 


By using this molar mass as well as the coefficients from Equation (B.5) and 
inserting these values into Equation (B.6) yields the specific capacity of NMC as 


ONMC, spec = —278mAhg !. (B.6) 


96.46 


B.2 Density from X-ray data 


The gravimetric density of a material i can be calculated by x-ray data. 
From [154], the density can be calculated by 


NM; 


= ————— B.7 
ЛА У; сеп E 


Oi 
where N is the number of atoms per unit cell, M; is the molar mass, Л/А is the 
Avogadro constant and V; сец is the unit cell volume. 


In the following, Equation (B.7) is applied to the positive electrode active 
material LiNi; /3Mnj /3C01/302, which is used in this work. From Appendix В.І, 
the molar mass M; = Mini Mni/5Coi/502 = Mwwc is computed as the 
stoichiometrically weighted sums of each element using their atomic weights. 
From Table В.І, the sum molar mass is 96.46g mol~!. Additionally, the values of 
the number of atoms per unit cell and the unit cell volume are taken from [141], 
where N = 3 and У; сеп = 100.769 A Note that these values refer to a fully 
lithiated NMC metal oxide, ie. y= 1 in LiyNi;/3Mnj/3Coj/50». However, 
densities for other lithiation states are computed accordingly. Finally, the 
gravimetric density of NMC is 


3.96.46 


=4770Ке m ?. B.8 
6.022 14076. 102 . 100.769 N Oe 


ONMC-— 
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